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ABSTRACT 


A mobile radio environment places fundamental limitations on the performance of 
wireless communication systems. Most models developed to predict propagation path 
loss have been historically performed in a statistical approach. These models are 
expensive to develop and do not offer the accuracy, computational advantages, and 
sufficiency as the parabolic equation (PE). The goal of this thesis is to develop a 3D 
model based on PE for predicting propagation path loss in urban areas on flat and hilly 
terrains. The PE method offers the computational advantages, where one can 
approximate the elliptic operator governing the true wave behavior by a much simpler 
parabolic operator that permits marching in range. Moreover those all-important aspects 
of propagation such as reflection and diffraction are included automatically in the 
formulation. Eour test problems on flat terrain and two test problems on hilly terrain will 
be simulated. Eor the flat terrain, the 3D PE model results will be compared with the 
two-ray, the four-ray, the UTD, and the numerical integration technique results. Eor the 
hilly terrain, the results of the 3D PE model will be compared with the UTD and the 
numerical integration technique results. 
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EXECUTIVE SUMMARY 


A mobile radio environment places fundamental limitations on the performance of 
wireless communication systems. Most models developed to predict the propagation path 
loss have been historically performed by statistical approach. These models are often 
complicated and expensive to develop and do not offer the accuracy, computational 
advantages, and sufficiency such as the parabolic equation (PE) technique or a ray based 
technique. 

The aim of this thesis is to develop a 3D model based on PE for predicting the 
propagation path loss in urban areas on flat and hilly terrains. The PE method offers the 
computational advantages, where one can approximate the elliptic operator governing the 
true wave behavior by a much simpler parabolic operator that permits marching in range. 
This method has the advantage that all-important aspects of propagation such as 
reflection and diffraction are included automatically in the formulation. Two types of 
terrains are considered, the flat earth and the hilly terrain. Eor the flat earth case, four test 
problems are examined. We compare the results from these cases with the results 
available in open literature from the two-ray model, the four-ray model, the UTD, the 
theoretical model proposed by Eee, and the numerical integration technique presented by 
Bertoni. Eor the hilly terrain case, two test problems are considered. The results of these 
cases are compared with the results presented by Piazzi and Bertoni. 

Eor the flat earth, the first test problem is to simulate and determine the 
propagation factor F at the final range over flat earth without obstacles being placed 
between the transmitter and the receiver. F is defined as the normalized field 

H^,[x,y,z)l where is the free-space magnetic field. A 

transmitting antenna height of Fit = 30 m, and the operating frequency of 1 GHz is used. 
The results of the 3D PE model are compared with the two-ray model. This is the 
simplest validation case we consider. 

The second test problem is to simulate and determine F with a single absorbing 
screen placed between the transmitter and the receiver that has a height Hk = 50 m and a 
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width, Wi: = 49.5 m. The screen represents a building. The transmitting antenna height is 
Ht = 60 m, and the operating frequency is 1 GHz. We compare the results of the 3D PE 
model with the results of the four-ray model which considers reflection and diffraction in 
the vertical plane joining the transmitter and receiver. 

In the third test problem, nine absorbing screens of uniform heights, equal 
spacing, variable screen widths are placed between the transmitter and the receiver. A 
transmitting antenna height of Ht= 10 m, and operating frequency of 1 GHz is used. This 
test case is equivalent to the radiowave propagating over a row of houses in an urban 
area. The three different screen widths of 50 m, 25 m, and 12.5 m are considered. Each 
screen has a height of 10 m and is spaced 50 m from the previous one. Their effects on 
the propagation factor F are studied. The propagation factors are determined at the 
rooftops. The results of the 3D PE model are compared with the UTD results presented 
by Andersen and the theoretical results proposed by Eee. 

The last test problem for the flat earth case involves 120 absorbing screens of 
uniform heights and equal spacing between the transmitter and the receiver. These 
screens represent row of buildings or houses. The transmitting antenna has a height of Fit 
= 125 m, and the operating frequency is 900 MHz. The screens have the height, Hk = 20 
m, and the width, Wk = 50 m. They are spaced 50 m apart. The propagation factor is 
determined at the final range. We compare the results of the 3D PE model with the 
results of the UTD method and the numerical integration technique presented by Bertoni. 

Eor the hilly terrain case, the first test problem involves a single rounded hill with 
multiple absorbing screens of uniform heights, equal spacing and variable screen widths. 
The three variable screen widths are 25 m, 50 m and 100 m. The screen height //yt = 7 m, 
and the spacing of 50 m are maintained constant. Their effects on the propagation factor 
F are studied. The transmitting antenna height is Ht = 57 m, and the operating frequency 
is 900 MHz. The field strengths are determined on the rooftops. We compare the results 
of the 3D PE model with the results of the numerical integration technique proposed by 
Piazzi and Bertoni. 
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For the second test problem in this class, the propagation took place over two hills 
of sinusoidal shape with multiple absorbing screens of uniform heights, equal spacing, 
and variable screen widths. The three variable screen widths are 25 m, 50 m and 100 m. 
A screen height = 7 m, and spacing of 50 m is used. Their effects on the propagation 
factor F are studied. The transmitting antenna has a height of Fit = 57 m, and the 
operating frequency is 900 MHz. Again, we compare the results of the 3D PE model 
with the results presented by Piazzi and Bertoni. 

Finally, the 3D model results based on PE for predicting propagation path loss in 
urban areas over flat and hilly terrains are simulated. Six different test problems are 
considered. We compare the results of the 3D PE model with the available results from 
the literature, and they show excellent agreement. We also demonstrate that the 3D PE 
model can support both flat and hilly terrains with multiple absorbing screens of uniform 
heights, equal spacing, and variable screen widths. 
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I. INTRODUCTION 


A. BACKGROUND 

A typical mobile radio environment in an urban area for wireless communication 
systems has no direct line-of-sight path between the transmitter and the receiver. The 
environment is so dynamic and the path between the transmitter and the receiver can vary 
drastically from simple line-of-sight to one that is severely obstructed by buildings, 
mountains, and trees. Hence, the mobile radio environment places fundamental 
limitations on the performance of wireless communication systems. Since the 
environment is extremely random, most analysis done on it has been historically 
performed in a statistical approach, often based on measurements made specifically for 
an intended communication system or spectrum allocation [Ref. 1]. 

Currently, there are large-scale propagation path loss models being utilized in 
commercial cellular communication systems in America, Europe, Asia, and around the 
world. Some of these popular models include the Okumura model, the Hata model, and 
the COST-231-Walfish-Ikegami model [Ref. 2]. As previously mentioned, these models 
were developed by measurements and statistical analysis made specifically for frequency 
coverage between 150 MHz and 2 GHz, and are typically based on a certain city 
environment like Tokyo or New York. These models provide reasonable approximation 
for current cellular communication applications. However, they are complicated and 
expensive to develop and do not offer the accuracy, computational advantages, and 
efficiency of models such as the parabolic equation (PE), ray methods, etc. 

Previously developed 2D parabolic equation based on the vertical plane method 

does not account for the lateral propagation of waves. Consequently, the mean path loss 

is overestimated and the standard deviation of error tends to be rather high [Ref. 3]. 

Therefore, a 3D PE formulation based on the split-step algorithm was developed and 

demonstrated for forward propagation around perfectly reflecting obstacles located on 

ground that includes the laterally propagating waves [Ref 4]. The 3D formulation is 

expected to substantially improve the accuracy relative to the 2D vertical plane method 

for the test problems. The goal of this thesis is to develop a 3D model based on a scalar 
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parabolic equation for predicting propagation path loss in urban areas on flat earth and 
hilly terrains. 

B. RESEARCH DEFINITION 

Given an operating frequency /, the transmitting antenna height Hi and the 3-dB 
beam width, and the terrain profile, we want to calculate the field at a certain range from 
the transmitting antenna. We use the parabolic equation (PE) method to determine the 
field at the desired location. Although there are several computational techniques for 
predicting the field at a desired range from the transmitting antenna including the ray 
tracing approach, and the integral equation method etc., few offer the computational 
advantages of the parabolic equation (PE) method, where one can approximate the 
elliptic operator governing the true wave behavior by a much simpler parabolic operator 
that permits marching in range. The PE method has the advantage that all-important 
aspects of propagation such as reflection and diffraction are included automatically in the 
formulation. However, the penalty for employing the PE method is that it neglects back 
scattering [Ref. 5]. This assumption will not contribute any significant errors for the 
class of applications considered in this thesis since radiowaves predominantly propagate 
in the forward direction. 

We consider a vertically polarized current source f{z) that produces only 

vertically polarized electric field component or equivalently a circumferential 
magnetic field component. We look at the y component of the magnetic field Hy because 
the depolarization is ignored. Thus initial magnetic field H{Q,y,z) at the location x = 0 
is set up by the current source on a flat plane, and it is directly related to the current 

source density J . It is easier to relate the magnetic field to the current source [Ref. 3, 4]. 
Eor this reason, we use the magnetic field in the 3D PE algorithm instead of electric field; 
however, the electric field and magnetic field in the far-zone are related via the medium 
impedance. The field is then split into even and odd parts to ensure that there are no 
discontinuities at zero crossing which leads to erroneous results during numerical 
evaluation. 
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In the presence of vertical screens that are perpendicular to the preferred axis {x- 
axis), we still ignore the depolarization. These screens represent a row of buildings or 
houses in an urban area. In this thesis, we assume that they are perfectly absorbing which 
means that all incident fields on them are zeroed out. Figure 1 illustrates how the 3D PE 
model algorithm handles the incident fields on the screens. First, we determine the fields 
H [x,y,z) in the absence of the screens and then make local adjustments to them. After 

making the adjustments, the field //_y,z)continues to march in range where the 

superscript indicates the fields immediately after the screen. This is then repeated till 
the desired receiver location is reached. 


Original Field 


Null Field 



> Y 


Figure 1. Fields Incident on a Perfectly Absorbing Screen. 

Therefore, this work will attempt to design a more accurate large-scale 
propagation path loss model based on the scalar 3D parabolic equation (PE) over the 
frequency band of 300 MHz to 10 GHz. The algorithm will predict the lateral and 
vertical wave propagation in an urban area that is comprised of vertical buildings with 
arbitrary cross section and perfectly absorbing surfaces on flat and hilly terrains. The 3D 
PE algorithm also includes the log-normal shadowing effects. We will compare the 3D 
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PE results to the two-ray model [Ref. 2], four-ray model [Ref. 2], the results given by 
Andersen [Ref. 6] obtained from the uniform theory of diffraction (UTD) method, the 
theoretical results given by Lee [Ref. 7], and the numerical results obtained by Piazzi and 
Bertoni [Refs. 8-10]. The 3D PE algorithm is implemented in MatLab. 

In Chapter II we provide detailed discussions on the numerical evaluation of the 
3D parabolic equation and the MatLab implementation. The current source has vertical 
extent along the z-axis and has a Gaussian distribution. The source standard deviation, 
(J^, is chosen to produce the required transmitting antenna 3 dB elevation beam width, 
0^;. A detailed discussion on the current source and its relationship to the magnetic field 
is also given in this chapter. 

In Chapter III and IV, we discuss our results and comparison efforts in detail. We 
compare the 3D parabolic equation (PE) results with the results reported in the literature. 
In Chapter IB, first, we examine a flat earth and perfectly reflecting ground. The results 
of the 3D PE model of this case are compared with the results of the two-ray model. 
Then we introduce single and multiple perfectly absorbing screens of uniform heights and 
equal spacing between the transmitter and the receiver over a flat earth and non-perfectly 
reflecting ground. These absorbing screens represent row of buildings of arbitrary 
widths. The effects of finite screen widths on simulation are also examined. 

In Chapter IV, we introduce a rounded hill and two hills of sinusoidal shape with 
multiple absorbing screens the uniform heights and equal spacing between the transmitter 
and the receiver. These scenarios are representative of well built-up urban areas. In both 
cases, we also examine the effects of finite screen widths versus the infinitely long screen 
widths on the simulation results. Results from these cases are compared with the results 
presented by Piazzi and Bertoni. Chapter V provides the conclusions and 
recommendations. All MatLab codes are listed in the Appendices. 
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II. NUMERICAL EVALUATION OF 3D PARABOLIC 

EQUATION 


A detailed method for evaluating the 3D parabolic equation (PE) numerically is 
presented in this chapter. First the 3D PE model algorithm flow diagram is presented in 
Section A. Then the initial current source and its characteristics are discussed in Section 
B. This is followed in Section C by a detailed discussion on how the magnetic field is 
related to the current source. A Hanning window is used throughout the algorithm to 
contain the fields both in spatial and wavenumber domains before taking 2D (Ny x Nz) 
Fourier transforms (IFFT’s) or 2D inverse Fourier transforms (IFFT’s). The Hanning 
window is discussed in detail in Section D. Brief descriptions of the auxiliary parameters 
are provided in Section E. 

In the parabolic equation the fields are coupled causally from one range to the 
next. There is no path for waves to return information from obstacles located ahead of 
the receiver. The initial fields are set up by the current source on an impedance plane. 
Presence of a vertical screen perpendicular to the preferred axis (v-axis) is handled by 
first determining the fields in the absence of the screen and then making local 
adjustments to them. This process is then repeated till the desired receiver location is 
reached. 

Figure 2 illustrates the geometry of the scalar 3D PE model approach. To find the 
field over an impedance plane at a particular range, x + Ax, given the field at a previous 
range, x, the latter is first decomposed into a spectrum of plane waves traveling in the 
positive v-direction. Each plane wave is then propagated to the new range using the 
appropriate propagator. The plane waves at the new range are then added to compose the 
field at X + Ak. The various spectral decompositions are carried out using 2D Fourier 
transforms (IFFTs) in the y-z plane. Furthermore, the magnetic field can be defined in 
terms of odd and even parts to ensure that there are no discontinuities at zero crossing 
which leads to erroneous results during the 2D Fourier transform (FFT) computation. In 
this thesis, the current source is assumed to be along the positive z-direction, and the 
antenna is vertically polarized. Over flat plane, the only non-zero field components are 
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E,,Ep, and Where p, z, the standard cylindrical coordinates, and ^is measured 
from the x-axis. We ignore depolarization and continue this assumption to hold in the 
presence of the screens. Thus the vertically polarized current source f {z) produces 

only vertically polarized electric field component or equivalently a circumferential 
magnetic field component. In this thesis, we choose to use the y component of the 
magnetic field Hy because the antenna depolarization is ignored. 



Figure 2. Geometry of the Scalar 3D PE Model [From: Ref. 11]. 

A. 3D PARABOLIC EQUATION ALGORITHM 

This section provides an overall view of the 3D PE model algorithm flow diagram 
as illustrated by Figure 3. The method of implementing the algorithm in MatLab is also 
presented. Obstacles’ placements and terrains are also discussed in the context of 
simulation. Additional details are provided in the subsequent sections. 
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Figure 3. 3D PE Numerical Evaluation Elow Diagram. 


The initial magnetic field H[0,y,z) is assumed to be generated by the current 


source on a flat plane. The Gaussian current source along the z-axis gives rise to a 
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circumferential magnetic field with a lateral component, H^.[x,y,z) ■ It is described in 

Section C. The transformed magnetic field H y {^0 ^) provides the initial field for 

the beginning of the main loop of the algorithm as shown in Figure 3, where the tilde 

indicates the transformed. When Hy[pi^,k^,k^ is multiplied by the free-space 

propagator , the initial transformed field has marched from the previous location x 
to a new location x + Ax and is given by equation (2.1). The field is split into even and 
odd parts for numerical convenience. 

H ^,{^x +Ax,ky,k^^ = < (2.1) 

The quantity Fn [k^) is the reflection coefficient of a plane wave for parallel polarization. 

Then the magnetic field H^ [x, y, z) at a new location x + Ax just right before the 
screen is found by double integrating equation (2.1) and is expressed by equation (2.2). 

(v,y,z) = -^ j Hy [x + Ax,k^,k^)e^''^"^''^^^dk^dk^ (2.2) 

—oo —oo 

If we examine equation (2.2), it is basically a 2D inverse Fourier transform, (IFFT), of 
the left hand side of equation (2.1) with respect to ky and k^. Then instead of evaluating 
equation (2.2) by complex double integrations, the 2D {Ny x NJ IFFT is used to evaluate 
it. The equivalent MatLab expression to equation (2.2) is defined by equation (2.3). 

M, f ' \ 

= (2.3) 

[27r) V J 

Also, MatLab divides the results by Ny x when it performs a 2D IFFT operation; 
therefore, the results are de-normalized by multiplying them by Ny x N^. The results are 
multiplied by a Hanning window in spatial domain to prevent aliasing effects. 
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At this point, after the magnetic field {x,y,z) has been found just before the 

screen, the obstacles are introduced that represent buildings, houses, trees, etc. in the 
simulation. In all the test problems, screens are used as obstacles and assumed to be 
perfectly absorbing. These screens represent rows of buildings having height Hk and 
width Wk. The Hk and 114 parameters are used to study the effects on the propagation 
factor and field strengths. All incident fields on the screens are eliminated since they are 
perfectly absorbing. The fields that propagate forward consist of the diffracted and the 
reflected waves. The H,y,z) indicates the field after adjusting for screens at 

location x. Then the fields begin their march immediately after the screen. Figure 4 
illustrates an example of a hilly terrain with multiple screens of uniform heights and 
equal spacing between the transmitter and the receiver. The propagation factor F and the 
field strengths may be determined at each marching step as the fields propagate forward. 
The data are stored as an array. 



Figure 4. Hilly Terrain and Multiple Equal Height and Equally Spaced Screens. 

Once the magnetic field H^{x,y,z), at a new location v, has been found, it may 

be expressed in terms of odd and even parts to ensure there are no discontinuities at zero 
crossing that cause erroneous results during the 2D Eourier transform (FFT) operations. 

The odd H (x, y, z) and even (x, y, z) fields are given as: 
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HAx,y,z) z>0 
-HAx,y,-z) z<0 


Hye[x,y,z) = 


HAx,y,z) z>0 
HAx,y,-z) z< 0 


The transformed quantities //>0 (j:, ) and Hye{^x,k^,kA are found by the 


following equations: 


HyQ[x,k^,kA= £° I” H ^.Q{x,y\z')e 
H ye[x,ky,kA = ^_^^_^H^A^,y\z')e 


Equations (2.6) and (2.7) are essentially the 2D Fourier transforms, 2D FFT’s, of 
Hyo (-t:, y,z) and//y^ (x, y, z) with respect to y and z. The equivalent Matlab operations 
of equations (2.6) and (2.7) are defined as: 

H yo[x,ky,kA = ^y^^FTli^H yQ{x,y,z)) (2.8) 

H ye (x,ky,kA = AyAzFFT 2(77^^ (x, y,z)) (2.9) 

Next the odd and even transformed magnetic fields are multiplied by the 
Hanning window in frequency domain to prevent aliasing. Then the odd transformed 

Hyo[x,ky,kAis multiplied by a factor—|^l-r||(A:^)^, and the even transformed 

FIye{^x,k^,,kA is multiplied by ^l^l + Fn . These fields are then combined to 
provide the newly transformed magnetic field immediately after the screen, i.e. 

H y [Q* ,ky,k^^. The steps in the main loop as illustrated in Figure 3 are repeated until 
the desired receiver location is reached. The propagation factor F at the final range is 
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defined as the normalized field H^.[x,y,z )//fy (/.5.)|, where is the free- 

space magnetic field and is defined as: 

1 ~ 

^v(/“^-) = -^’^oSin(^)cos(0)—^(^ocos(0))-— (2.10) 

4-71 K 

where 0and ^are the angles in spherical coordinates, i = ^I^, and g{k^) defined by 

equation (2.20) is the vertical plane response of the current source, is a range between 
the transmitter and receiver, and ko is the free space wave number. R, p, 9, and ^ are 
defined by the following equations: 



(2.11) 

p = ^x^ + y^ 

(2.12) 


(2.13) 


^ = sin — 

(2.14) 

[PJ 



The current density and its characteristics are discussed in the next section. 

B. GAUSSIAN CURRENT SOURCE 

The fields are excited by a vertically polarized current source with a prescribed 
aperture distribution / (z). We choose a Gaussian distribution current source to allow us 

to match the source standard deviation, O ^, to the transmitting antenna 3 dB elevation 
beam width, 0^^ [Ref. 3, 4]. It is also convenient to relate the magnetic field to the 
current source. We consider a transmitting antenna at height z = Hi and the current 

source density J as defined by equation (2.15). 

J = zIjS{y-0)f{z) (2.15) 

where 
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1 


(2.16) 


/(^) = - 


Az-H,fl 2<J^ 




is a Gaussian current source with a standard deviation of (J^, IJ. is a current moment 
[Ref. 11], and —O) is a dirac function. Figure 5 shows the vertically polarized 

Gaussian current density function J placed over an impedance plane. The current is 
assumed to be independent of y and has a peak at z = //,. is its standard deviation and 
may be chosen from Table 1 to fit the required transmitting antenna 3 dB elevation beam 


width, 0^; [Ref. 3]. It may be observed that R = — + /?^ is the direct distance 

V 2 2 

X + y is the horizontal distance from the source, and 0 and (j) 


are angles in spherical coordinate defined by ^ = sin ^ — and (j) = sin ^ 


^ P? 

\Rj 


P) 


. For the 


geometry shown, the only non-zero field components are , Ep, and H, 


/I 



Figure 5. Gaussian Density Current Source over an Impedance Plane. 
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Figure 6 illustrates the Gaussian current source function/fzj with Ht= 10 m and 


(J^ = X. We discuss the magnetic field in detail in the next section. 


®ei ’ deg 

(7^ IX 

60 

0.20 

45 

0.30 

30 

0.45 

15 

1.00 

10 

1.52 

5 

3.03 


Table 1. Transmitting Antenna Elevation Beam Width Versus t7 [From: Ref. 31. 










The function / (z) can be written in terms of its odd and even functions, (z) 
and/g (z), respectively, to permit a more accurate computation of the Fourier transform 
in the vertical direction. The odd and even extensions of /(z) are given in equations 


(2.17) and (2.18). 



j/(z) z>0 

|-/(-z) z<0 


fe{z) 


j/(z) z>0 

|/(-z) z<0 


(2.17) 

(2.18) 


The functions/^ (z) and /^(z) are shown in Figures 7 and 8. If we add the/^(z) 

and/g(z) together and divide by Vi, we will get back the original function f{z). 

Hence, the function decomposition into its odd and even parts to facilitate computation is 
a valid one. 



fo(^) 

Figure 7. Odd Part of / (z). 
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Figure 8. Even Part of f {z). 


We need to examine the Fourier transform of the current source density J . We 
shall use the transformed information to derive the initial transformed magnetic field in 

Section B. Since J is independent of y, if we take the Fourier transform of equation 
(2.15) with respect to z, which is tantamount to taking the Fourier transform (FFT) of 
f {z) ■ The FFT of J(_y —O) with respect to y is equal to 1. Therefore, the Fourier 

transform of f{z) with respect to z is given below: 


fiK) 


-ikH, - 

e ^ 'e 


s{K)e 


(2.19) 
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It is easy to see that 


fe{K) = ^s{K)^o^{KHz) 

(2.21) 

fo{k,) = -2ig{k,)^m{k^H,) 

(2.22) 


C. MAGNETIC FIELDS 


In this section, we derive the initial transformed magnetic field Hy{^x,k^,k^ 


from the transformed current source f {k^). The fields march in the transformed space, 

and the propagation factor or the field strength is determined in the physical space. Then 
the basic 3D parabolic equation (PE) algorithm will start with the initial transformed 

magnetic field H y[Qi^ ,ky,k^^ as illustrated in Figure 3, where 0“^ indicates the field 
immediately after the screen. Since there is no screen for the initial field, 0"^ is equal to 0. 
First we need to define the relationship between the initial magnetic field 

H^ y,z) and the current density function J . The initial H^ y,z) relates to the 


current density J is described by equation (2.23) [Ref. 11] 

^ -oo 0 

- f [r, (<:,)] ] dy ~\dz 'S(y0)f(z y‘~’' 

“ (2.23) 

T ] (t )d*, -;2;Res[r. (i ,)] 
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where the second part of equation (2.23) describes the surface waves. In many long- 
range propagation problems, particularly in the VHP band and above, it is possible to 
completely ignore the surface wave terms as they decay rapidly with range [Ref. 5]. 

Then equation (2.23) is reduced to a simple relationship between the H^ (0’^,_y,z) and 

the Gaussian current source function f [z]. 

Again, to ensure there is no discontinuity at zero crossing that gives erroneous 

results in the 3D PE numerical calculation, the magnetic field may be defined in terms of 
its odd and even parts, [x,y,z) and H (x, y, z ), respectively. The methods are the 
same as in equations (2.4) and (2.5). 

From equations (2.4), (2.5), (2.23), and the transformed odd fo{k^) and even 
f current sources as defined by equations (2.21) and (2.22), we can easily derive 
the initial transformed odd Hyo{^0 ^and even Hye{pi^,ky,k^^ magnetic fields. 
The initial odd transformed magnetic field H yo , k^ ,k^^ is defined by equation (2.24) 

Hyn(0\k,,k^)=^fJL) 

\ y z; 2 (2.24) 


The initial even transformed magnetic field Hye{^0^,ky,k^^ is defined by equation 
(2.25) 

V y z/ 2 (2.25) 


These fields are defined as column vectors having elements. From Section A, we 
know that the current density function J is independent of y, and the 
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FFr(j();-0)) = l; then we can repeat the column of the initial 

transformed H > 0 (O"^,) and Hye{pi^,k^,k^^ Ny times. Now we have the 2D initial 

transformed/f>o(0'^,A:^,A:^) mAHye{Qi^,ky,k^^ with dimension of Ny x N^. The 
number of elements along the y-axis and z-axis are Ny and N^, respectively. 


The transformed magnetic field Hy {^x,k^,,k ^) is defined by equation (2.26). 

Hy(x,k^,k^)= (2.26) 

0 


Recall that e and are the free-space propagators for direct and reflected rays, 
respectively, and y’ and z’ are dummy variables of integration. r||(A:J is the reflection 
coefficient for the parallel polarization as a function of kz in the complex kz plane and in 
terms of normalized impedance Zs. r,, [k ^) and Zs are expressed as: 


r„(^J 


K-K^s 

K + K^s 


= 



(2.27) 

(2.28) 


e =e + 


)18(j(m5 / m) 
f(MHz) 


(2.29) 


where/is the operating frequency measured in MHz. We use equations (2.4) and (2.5) to 
decompose the magnetic field H^.[x,y\z') into Hy^[x,y\z') and H^,^[x,y\z') and 

substitute them into equation (2.26). Now we add theH^^ (x,y',z')and H^.^{x,y',z') 

together for z > 0 and perform variable substitutions to change the inner limit of 
integration from 0 to 0 ° to - 0 ° to 0 ° and divide it by ¥ 2 , and similarly for z < 0. The result 
is: 
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H,{x.k,Xy\[l + T,(kS\[_[^H_„(x.y\z')e 


-i[k y'+k^z'] 


dy'dz' 


_1 -r,,J] £ £^^, 0 [x,y\z')e 


-i(k^y'+k^z'} 


dy'dz' 


(2.30) 


= ^ [i+r,, (^ J] -[i - r„ (J]//,0 ( j:,J 


If we examine equation (2.30), we see that the double integrations of 
Hye{x,y',z') and H^^{x,y',z') are the 2D Fourier transforms (FFT’s) of these fields. 
Thus equation (2.30) can be evaluated with the 2D (Ny x Nd FFT’s instead of complex 

double integrations. Thus H>• £,)can be defined in terms of FIyo[x,y,z) and 


Hye{x,y,z). This is how we are going use the 2D FFT’s in the MatLab codes to 
implement equation (2.30). 

After making the adjustment for the screens, the fields begin their march in range 
immediately after the screen expressed as 


Hy[x^,k^,k^ 




(231) 


Equation (2.31) provides the initial transformed magnetic field for the beginning of the 
main loop of the 3D PE algorithm as shown in Eigure 3. Once again, indicates the 

fields immediately after the screen. The initial fields FIyo{^x^,k^,,k^^ and 

H yei^x^ ,ky,k^^ are defined as matrices having Ny x elements. These fields will allow 

the algorithm to account for all the vertical and lateral propagation of waves, which is the 
fundamental goal of the 3D PE formulation. Previously developed 2D PE algorithms 
based on the vertical plane method did not account for the laterally propagating of waves. 
Consequently, the mean path loss was overestimated and the standard deviation of error 
tends to be higher [Ref. 11]. Therefore, the 3D PE formulation is expected to 


19 



substantially improve the accuracy relative to the 2D vertical plane method for the test 
problems. 

When we previously performed column-repeating the initial H 
andH,k^,k^^ by Ny times to create the 2D (Ny x N^) H yo{Qi^ ,ky,k^^ and 2D 

Hye{Ql^,k^,,k^^, it was to ensure that the matrices always have dimension of Ny x 

since this is a general aperture case. When Ny is not equal to N^., we have a rectangular 
aperture, and when they are equal to each other, we have a square aperture. The 3D PE 
algorithm is implemented to support any Ny x Nz aperture. Where Ny and Nz are elements 
along the vertical and horizontal ranges, respectively, with a suitable power of 2, i.e., Ny 
= 1024 or 2048. 

Discussion thus far pertains to what happens outside of the main loop of the 
algorithm. Equation (2.31) provides the initial transformed field for the beginning of the 

main loop for the 3D PE algorithm; then we multiply H ,k^,k^^ by the Hanning 

window in frequency domain before marching. We always multiply the fields by the 
Hanning window before taking 2D Eourier transform (FFT’s) or 2D inverse Eourier 
transform (IFFT’s) to contain them in space and wavenumber domains. The Hanning 
window is discussed in the next section. 

D. THE HANNING WINDOW 

In order to contain the fields in space and wavenumber domains we multiply them 
with the Hanning window [Ref. 3] before taking 2D Eourier transform (FFT’s) or 2D 
inverse Eourier transform (IFFT’s). The Hanning window provides a gradual rolloff to 
zero over the last quarter of the domain, and may be constructed from the two Hanning 
sequences below. Equations (2.32) and (2.33) are the Hanning sequences along the 
horizontal range (y-direction) and the vertical range (z-direction), respectively. The 
mirror images of and h^{t^ below t = 0 are used for negative wave numbers. 
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for 0<?<3AA^,/8 


n 


K{t) 


1 sin^ 




V ^.v y 


for 3NJS<t<Nyl2 


Kit) 


sin 






for 0<?<3A^^/8 
for 3NJ%<t<NJ2. 


(2.32) 


(2.33) 


The Hanning window, in general, has a rectangular shape since Ny and are not 
necessarily equal. When Ny and are equal, then the Hanning window has a square 
shape as discussed in the previous section. Figures 9 and 10 illustrate the examples of 
512 elements of \ [t) and [t) Hanning sequences. 



Figure 9. K {t) Hanning Sequence. 


21 






















500 



Figure 10. [t) Hanning Sequence. 

As in the previous section, the Hanning window must have a dimension of Ny x 
N^. Again, the procedure is similar to Section C; we first create a matrix from the 
Hanning sequence h^{t) by column repeating the sequence Ny times. This method 

creates a. Ny x matrix from the h^{t) sequence. Then we construct the second matrix 
from the Hanning sequence the \{t) by row-repeating it N^ times (thereby constructing 
a Ny X Nz matrix from h^.{t) sequence). We column-repeated \{t) Ny times to also 
ensure that we do not create a matrix that has a column dimension greater than Ny. The 
same is true for the h^, (f) case. We row-repeated N^ times to ensure that the row of the 

matrix does not have dimension greater than Nz- This implements the Hanning window 
in MatLab. 

Figure 11 shows the two Hanning matrices created from (f) and h^. (f). To 

construct the final Hanning window, the two matrices are multiplied together element by 

element. Figure 12 shows an example of a 512 x 512 Hanning window. A similar 
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Hanning window is used in the 3D PE algorithm depending on the dimensions of Ny and 

Nz. 



Figure 11. Ki^) (f) Hanning Matrices. 
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Figure 12. 3D Hanning Window. 
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E. SPATIAL AND FREQUENCY PARAMETERS 

The values of Ak, Ay, Az, Aky^ and Ak^ and the range of kx, ky, and k^ are dictated 
by the Nyquist sampling theorem. Ax is the range increment, Ay is the horizontal 
increment, and Az is the vertical increment, kx , ky, and k^ are the range wavenumber, the 
horizontal wavenumber, and the vertical height wavenumber, respectively. Aky and Ak^ 
are the increments of the horizontal and the vertical wavenumbers, respectively. 
References 3,11, and 12 provide details to obtain the values for these parameters. 

For all the test problems in this thesis, we assume that Ay = Az = X and Oz = X. 
We pick Ax between 25 m and 150 m. If we choose a smaller Ax, the computation takes 
longer. If we choose a larger Ax, the computation is accelerated, but we might miss the 
obstacles that are within in the marching step. Therefore, Ax has to be optimally chosen 
accordingly for each problem. However, Ax can be chosen to have variable values within 
a simulation. Furthermore, for Oz = X corresponding to the transmitting antenna of a 3- 
dB bandwidth of 15° is shown in Table 2. The value of Oz can be picked from Table 2 to 
match the transmitting antenna 3-dB bandwidth. 

We compute the propagation factor F at the final range as well as at each 
marching step. We also consider the field strengths at each marching step and on the 
rooftops. In Chapter HI, we shall present the results for the flat earth case. The hilly 
terrain test problems are presented in Chapter IV. 
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III. FLAT EARTH RESULTS 


In this chapter, we present results for the propagation factors for the flat earth and 
perfectly reflecting ground. The propagation factor at the final range without obstacles 
between the transmitter and the receiver is simulated in Section A. The result of the 3D 
PE model is compared with the result of the two-ray model [Ref 2]. Next, in Section B, 
we place a single absorbing screen between the transmitter and the receiver and present 
the result for the propagation factor at the final range. The screen represents a building. 
The simulation result is compared with the result of the four-ray model [Ref. 2]. Then, 
in Section C, we place nine absorbing screens of uniform heights, equal spacing, and 
variable widths between the transmitter and the receiver. These screens represent a row 
of buildings or houses in residential areas of a city. In this case, the relative propagation 
factors are measured at the rooftops for each marching step. The results of the 3D PE 
model are compared with the results presented by Andersen [Ref. 6] and Eee [Ref. 7]. 
Einally in Section D, we examine 120 absorbing screens of uniform heights and equal 
spacing and measure the relative propagation factor at the final range. Again, these 
screens represent a row of buildings or houses in residential areas of a city. The result of 
the 3D PE model is compared with the numerical integration technique result proposed 
by Bertoni [Refs. 10]. 

A. 3D PE AND THE TWO-RAY MODEL RESULTS COMPARISON 


In this section we stimulate the two-ray model and the 3D PE model with the 
same parameters and compare their results. Eigure 13 shows the two-ray model over flat 
earth. The two-ray model propagation factor, F 2 ray, [Ref. 2] is defined as follow with 


2 ray 


l + Te 


ihgAR 


(3.1) 


where A/? =-, T = , if/= idin 


d 


sin + 


^ 18(7(m5/m) 

, — £ -\-- 




\ d j 


fiMhz) 


and Z, = ^£^^-cos\y/)/£^^ - » 


0 . 
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The H, and Hr are the heights of the transmitter and receiver’s antennas, 
respectively, and d is the horizontal distance between the transmitter and receiver. A/? is 
the path difference between direct and reflected rays, and /"is the reflection coefficient 
for parallel polarization as a function of the grazing angle \|/ and the complex impedance 
Zs. and S are the complex dielectric constant and conductivity of the medium, 
respectively. It has been assumed that J » 7/^ and . 



Figure 13. Two-ray Model Over Flat Earth. 

Figure 14 shows the vertical cut of the relative propagation factors, F’s, of the 
two-ray model and 3D PE model. The two models are stimulated by the same parameters 
as listed on Eigure 14. The results were taken at a range of 1000 m from the transmitter a 
over flat plane. The circle with the solid line corresponds to the two-ray model and the 
plain solid line is the 3D PE result. The difference between the two results is within 1%. 
If we increase the height of the receiving antenna, the reflected ray will tend to decrease 
and the direct ray will dominate. Therefore, we expect the two F’s to approach zero dB 
for large heights. 
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Figure 14. 3D PE and Two-ray Comparison (Vertical Cut). 

B. 3D PE AND FOUR-RAY RESULTS COMPARISON 

In Section A, we have shown that the result of the 3D PE agreed with the two-ray 
results. Now, we place a 50 m by 49.5 m absorbing screen (a single knife-edge) in the 
3D PE model located 125 m from the transmitter and 375 m to the receiver with the 
frequency of operation of 1 GHz. The loss in this case depends on the height of the 
knife-edge above ground, its relative location from the transmitter/receiver, the ground 
constants, and the frequency of operation. Eigure 15 shows an absorbing knife-edge 
between the transmitter and the receiver. Table 2 provides the quantities required for the 
evaluation of the propagation factor F^. ^ for the 2D-four-ray model [Ref. 2]. The 3D PE 

model and the four-ray model are simulated for the same assumed parameters. The four 
paths may be identified from the transmitter to receiver via the tip of the knife-edge. The 
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total received signal is the sum of the direct wave, if any, and the diffracted waves 
received via each one of the four paths. The total diffracted field in the presence of the 
knife-edge is defined by equation (3.2). 

(3.2) 

n=\ 



Figure 15. An Absorbing Knife-Edge Between Transmitter and Receiver. 


Case n 

Transmitter 

Location 

Receiver 

Location 

Free-Space Path 
Length rn 

Clearance 
Height hn 

Reflection 
Coefficient, r„ 

1 

T, 

R. 


d 

1 

2 

T/ 

Rx 

TR^=^F+{H,+H^y 

d 

Fa 

3 

T, 

Rx’ 

tX = tX 

* d 

Fb 

4 

T/ 

Rx 

tX=tr. 

d 

FaFb 


Table 2. Quantities For The Evaluation of [After: Ref. 2] 
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where Eon is the free-space field for the situation, Fn is the normalized field for that 
path, and Hk is the height of the knife-edge. Assuming omni-directional patterns for the 
transmitting and receiving antennas the free-space field is defined as follows 

Eon= -,n=l,...,4, (3.3) 


where r„ is the free-space path length for the situation, and the individual knife-edge 
normalized fields are given by 



(1+j)^ 




r„, 


(3.4) 


where F is the Fresnel integral, Fn is the total reflection factor, and h„ is the clearance 


height for the path. Hi is the radius of the first ellipsoid in the plane perpendicular to 
the line-of-sight (LOS) path and can be obtained from equation (3.5): 


^1 = 


I 

r/j -f d2 


(3.5) 


where d; is the distance between the transmitter to the knife-edge, d 2 is the distance 
between the knife-edge and the receiver, and X is the operating wavelength. The 
propagation factor Fj. ^ of the four-ray model can then be obtained from equation (3.6). 


k.e. 



(3.6) 


Additional details of the four-ray model are provided in reference 2. 

Figure 16 shows the vertical cut of the relative propagation factors, F’s, of the 
four-ray model and 3D PE model. The circle with the solid line corresponds to the four- 
ray model and the plain solid line is the 3D PE result. The results were taken at a range 
of 500 m from the transmitter over flat plane. The screen was placed 125 m from the 
transmitter and 375 m to the receiver. The two results are in excellent agreement. As 
previously mentioned, the relative propagation factors F’s approach zero dB as we 
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increase the receiving antenna height. In this case, it is already close to zero dB at a 
receiving antenna height of only 80 meters. 



Figure 16. 3D PE and Four-ray Comparison (Vertical Cut). 

C. 3D PE MODEL AND ANDERSEN’S AND LEE’S RESULTS 
COMPARISON 

In the previous section, we have demonstrated that the results of 3D PE model 
corresponded well with the results of the four-ray model. In this section, we consider 
nine absorbing screens (multiple knife-edges) of uniform heights, equal spacing, and 
variable widths. The results of the 3D PE model are compared with the results presented 
by Andersen [Ref. 6] and Fee [Ref. 7]. Andersen’s results were obtained from the 
uniform theory of diffraction (UTD) method. Fee presented the theoretical results in 
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reference 7. Both the UTD and theoretical results assumed that the screen widths are 
infinitely long, and the screen heights are finite. Full details on the UTD method are 
provided by references 6, 7, and 14. 

We also examine the effects of three finite screen widths on the relative 
propagation factors, F’s. The three finite screen widths are 50 m, 25 m, and 12.5 m. In 
each case, nine absorbing screens were placed between the transmitter and the receiver 
with each screen having a height of //i:=10 m. The screens are placed 100 m apart from 
the each others. The first screen is 100 m from the transmitting antenna, and the last 
screen is 100 m to the receiving antenna. The distance between the transmitter and the 
receiver is 1000 m. The propagation factors are determined at the top of the screens. 


The theoretical propagation factor, Fth, in dB is defined by equation (3.7). It is a 
function of edge numbers, where the edge number is N +1, and N is a number of screens. 

^ I ^ 


20 log 


10 


A^ + I 


(3.7) 


Again, these screens represent a row of buildings or houses in residential areas of 
a city. In these types of environments, the base station antennas of cellular 
communication systems are typically located above or near to the rooftops of the 
surrounding buildings. In these cases, the propagation takes place over the buildings, 
which can be modeled by multiple forward diffraction past rows of buildings. We model 
each row of buildings as an absorbing knife-edge, and via the 3D PE numerical technique 
we determine the loss associated with multiple forward diffraction over the knife-edges. 
Furthermore, depending on their construction, buildings may have a flat roof, a peaked 
roof, a flat roof with a parapet, or a myriad of other roof designs [Ref. 13]. Figure 17 
illustrates building profiles that mobile communication systems designers can be 
expected to encounter for typical urban environments between the base stations and 
mobile units. Figure 17 (a) shows that two absorbing screens model a double parapets 
roof building. Figure 17 (b) illustrates that a single parapet roof building is modeled by 
one absorbing screen with the screen placement at the start of the building. And Figure 17 
(c) shows that a peaked roof building is also modeled by one screen, but the screen 
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placement is at the center of the building corresponding to the peak of the building. 
Reference 13 provides more details on building profiles and screen placements. 



Figure 17. Typical Building Profiles in Urban Areas, with Their Equivalent Screens 

Placements [After: Ref. 13]. 

Figure 18 illustrates the multiple absorbing knife-edges with equal heights and 
equal spacing for this test problem. The assumed parameters used in the 3D PE model 
are listed on Figure 19. 
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Relative Prop Factor (dB) j,;-- 



igure 18. 


Perfectly Absorbing Equal Heights and Equal Spacing Multiple Screens. 



Eigure 19. 3D PE, UTD, and Theoretical Results. 
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Figure 19 shows the results of the 3D PE model with the results presented by 
Andersen and Lee. The solid line is the theoretical result for 2D. The circle with solid 
line represents the results presented by Andersen. The triangle with dashed line are the 
results of 50 m screen width, the diamond with solid line are the results of the 25 m 
screen width, and the star with dashed line are the results of the 12.5 m screen width. The 
UTD, the 50 m and 25 m screen width cases, the curve slope track the theoretical closely 
up to the seventh screen. After the seventh screen, they are only 0.5 dB different from 
the theoretical. For the case of 12.5 m screen width, the results do not quite agree with 
the theoretical, UTD, 50 m, and 25 m results, but the difference is still less than 1 dB. 
These differences are caused by the fact that the smaller screen width allows the lateral 
waves and diffracted waves to constructively add or destructively add. These lateral and 
diffracted waves might have added constructively when they arrive at the receiving 
antenna up to the seventh screen; then added destructively after the seventh screen, which 
resulted in a 1 dB different from the theoretical and approximately 1.5 dB different from 
the UTD and the two cases of the 3D PE model. However, we expect to see the results of 
the 3D PE model approach the theoretical results if we consider wider screen widths that 
allow less lateral waves to contribute to the overall results. 

All presented data used absorbing screens of uniform heights and equal spacing, 
but the 3D PE algorithm is capable of supporting multiple screens of variable heights, 
widths, and spacing. Based on the comparison results in sections A, B, and this section, 
we may assume that the 3D PE model is capable of computing the propagation factors 
over flat earth with multiple absorbing screens of non-uniform heights, unequal spacing, 
and variable screen widths. Next we consider 120 screens of uniform heights and equal 
spacing over flat earth and perfectly reflecting ground. 

D. 3D PE MODEL AND 120 UNIFORM HEIGTH AND EQUAL SPACING 

SCREENS 

In this section we consider 120 screens of uniform heights and equal spacing. 
Again, the propagation takes place over the buildings. As previously mentioned, the 3D 
PE model is capable of supporting multiple absorbing screens of variable heights, widths, 
and spacing. But for the purpose of this study, we will only consider absorbing screens 
of uniform heights and equal spacing. We compare the 3D PE model results with the 
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results presented by Bertoni [Ref. 10]. Figure 20 shows setup geometry of this test case. 
We consider a plane wave incident at an angle a = 1° at an operating frequency of 900 
MHz. The distance between screens is 50 m, which means that the distance between the 
transmitter and the receiver is 6050 m apart. Based on the incident angle a of 1° and the 
separation distance between the transmitter and the last screen, we determine that the 
vertical clearance height, Hch, from the top of the 120* screen to the tip of the 
transmitting antenna is 105 m. If the screen height of 20 m is chosen, then the 
transmitting antenna height should be 125 m. We assume the screens width to be 50 m. 



Figure 20. 120 Equal Height and Equal Spacing Screens. 

Eigure 21 shows the height variation of the relative field strength, F, computed by 
3D PE model, incident on the row of 120 screens of 20 m heights and 50 m spacing for 
the parameters assumed the previous paragraph. The receiving antenna height is 
measured in wavelength (m) while the magnitude of the field strength is in linear scale. 
Eigure 19 also shows the simple diffraction in the shadow region of the receiving antenna 
Hr < 0. We use the rooftop as the reference point for Hr = 0. Above the rooftops, the 
field variation is similar to that of a standing wave resulting from the summation of the 
direct waves and reflected waves from the plane of the rooftops. The reflected waves are 
in fact the sum of the waves diffracted from the rooftops, whose phase variations cause 
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them to add constructively in the specular direction and other grating lobe directions. 
Because the diffraction coefficients decrease with angle, the reflected field is greatest in 
the specular direction. 


Ny = 1024 Nz = 1024 Dx = 25 meters Ht = 125 meters 



Receiving Antenna in Wavelegth 


Figure 21. 120 3D PE Results For Equal Height and Equal Spacing. 

Figure 22 shows the height dependence of the field strength, F, of the 3D PE 

model with the result from the numerical integration method presented by Bertoni and the 

result from the UTD [Ref. 10]. The circle with solid line is the sampling result of the 3D 

PE model shown in Figure 21. The asterisk with solid line is the result of the numerical 

integration method, and the triangle with solid line is the result of the UTD. In the 

shadow region, the results of the 3D PE model and the numerical integration method 

have good agreement, but the result of the UTD is slight higher, which implies that the 

UTD overestimates the field amplitude [Ref. 10]. Above the shadow region, the 
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IV. HILLY TERRAINS RESULTS 


In this chapter, the propagation takes plaee over hilly terrains, and we study the 
field strengths. A single round is eonsidered in Seetion A. The two hills of sinusoidal 
shape are eonsidered in Seetion B. In both oases, we place multiple absorbing soreens of 
uniform heights, equal spaeing, and variable widths between the transmitter and the 
reoeiver. We measure the field strengths, instead of propagation factors, at the rooftops. 
The results of the 3D PE model are oompared with the results presented by Piazzi and 
Bertoni [Refs. 8 - 10]. 

A. A ROUND HILL AND 3D PE MODEL 

In the previous test oases, we have examined the propagation faotors over flat 
earth with multiple soreens of uniform heights, equal spaeing, and variable widths. Now, 
we oonsider a rounded hill with multiple absorbing soreens of uniform heights, equal 
spaeing, and variable widths. We examine the diffraeted field strengths at the rooftops, 
instead of the propagation faotors as in the previous oases. We also oonsider the effeots 
of finite soreen widths on the field strengths. We oompare simulation results of the 3D 
PE model with results presented by Piazzi [Ref. 8]. Eigure 23 illustrates setup geometry 
for this test ease. These sereens represent a row of buildings or houses in residential 
areas of a well built-up eity. The buildings are 7 m high and the row separation Ds is 50 
m, and the soreen widths are Wt = 100 m, 50 m, and 25 m. The oylindrioal hill has a 
radius Rh of 10 km [Ref. 8]. The maximum height of the hill is 50 m. The base of the 
hill is at a distanoe of Xh =1000 m from the peak. Eor the simulation, the transmitting 
antenna is assumed to be loeated at the point x = -1000 m from the hilltop and has a 
height of Ht = 57 m, whieh plaees it at the same height as the highest rooftops [Ref. 10]. 
The eurrent source is exeited at an operating frequency of 900 MHz. The rounded hill is 
defined by equation (4.1). 

( 4 , 1 ) 

where Rh is the radius of the hill, Hk is the height of the buildings, and x is the range in 
meter. The total simulation range is 6000 m. 
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Figure 23. A Round Hill with Multiple Sereens of Uniform Heights and Equal 

Distanees. 

Figure 24 shows the profile of multiple absorbing screens representing Figure 23 
used to carry out the simulation in the 3D PE model. Furthermore, if we can model an 
urban area with multiple absorbing screens of variable heights, widths, and spacing, the 
3D PE model can evaluate the propagation factors as well as field strengths at any point 
in the simulation. 

Figure 25 shows the results of 3D PE model, the free-space result, and Piazzi’s 
result. There are three 3D PE model results. The solid line is the result of the 3D PE 
model with 100 m screen width. The dashed line is the result of the 50 m screen width. 
The dashed and dot line is the result of the 25 m screen width. The triangle with solid 
line is the free-space result, and the circles are Piazzi’s results. Up to 1000 m away from 
the base station, the field strength behaves like free-space. Then at the hilltop the field 
behaves like a single knife-edge case, where the magnitude of the field strength is 6 dB 
below the free-space. Between the base station and the peak of the hill, the field strength 
does not behave like a free-space field. Perhaps this is one of the limitations for the 3D 
PE model, but this is minor issue. After the hilltop and behind the hill up to x = 1000 m, 
the field strength decays rapidly and linearly with distance at a rate of approximately 
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0.055 dB per meter. This may be understood in terms of the creeping wave 
representation of the fields traveling over the houses. At the point of v = 1000 m from 
the hilltop, the field strength reaches its lowest point, but begins to increase rapidly to 
maximum and then decreases slowly with distance. We expect the field to continue 
decaying as it propagates away from the base of the hill, but that is not the case because 
the diffracted waves from the rooftops on the hillside. We also expect the field to 
eventually level off and begin to decay and behave like multiple knife-edges propagation 
case. As shown in Figure 25, the field levels off around x = 3000 m away from the base 
of the hill. If we continue to measure the field strength for both the 3D PE model and the 
free-space, the field strength of the 3D PE model will eventually be 6 dB below the free- 
space value as in the a single knife-edge case. 



Eigure 24. Side View of a Single Hill with Rows of Uniform Heights and Spacing 


Buildings. 
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Figure 25 also shows that we have excellent agreement between Piazzi’s result 
and the result of the 3D PE model especially with screen width Wit =100 m. However, it 
is not the case for the 50 m and 25 m screen widths. Furthermore, for both the 50 m and 
25 m screen widths, there are oscillations as the fields move away from the base of the 
hill. There appears to be more oscillation in the 25 m screen width case than for the 50 m 
case. Its magnitude is also higher than the other cases, but it begins to approach other 
cases as it moves away from the base of the hill on backside. The oscillations and higher 
magnitudes are caused by lateral waves because of the finite screen widths. Therefore, if 
we increase the screen width, then we expect the result 3D PE model to agree with 
Piazzi’s result because Piazzi assumed the screen width is infinitely long. This is 
equivalent to the 2D PE model which did not take into account the lateral waves. This is 
demonstrated the 100 m screen width case, which prevented the lateral waves from 
significantly contributing to the overall result. We have shown that the 3D PE model is 
capable of supporting the variable screen widths. The last test problem we will present in 
this thesis is the two sinusoidal hills. 



Screen Position (m) 


Figure 25. A Single Hill with Rows of Buildings of Uniform Heights and Spacing 

Field Strength. 
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B. 3D PE MODEL AND TWO HILLS OF SINUSOIDAL SHAPE 


In the previous section, we have considered a single rounded hill and 
demonstrated that the result of the 3D PE model agreed with Piazzi’s result. We have 
also shown that the 3D PE model can support variable screen widths. Now we consider 
two sinusoidal hills, and this is the last test problem we will consider for this thesis. It 
has similar setup geometry as shown in Eigure 23, except it has two hills of sinusoidal 
shape. Eigure 26 illustrates setup geometry for this test case. 



Eigure 26. Two Hills of Sinusoidal Shape. 


Again, these screens represent a row of buildings or houses in residential areas of 
a well built-up city. The buildings are 7 m high and the row separation Ds is 50 m with 
building widths Wk = 100 m, 50 m, and 25 m. The maximum heights of the two hills are 
50 m. The two sinusoidal shape hills have period of 3000 m [Ref. 9], which means the 
distance from peak to peak is 3000 m. The distance Xh from the first peak to its base, 
measured toward the second peak, is 1500 m. Eor the simulation, the transmitting 
antenna is assumed to be located at the point x = -1500 m from the first hilltop and has a 
height of H, = 57 m, which again places it at the same height as the highest rooftops [Ref. 
10]. The two hills of sinusoidal shape are defined by 


Z{x) = Hi^ -I-25-I-25cos 


^ 27rx ^ 
,3000, 


(4.2) 
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where Hk is the height of the buildings and x is the range in meters. 

Figure 27 shows the profile of multiple absorbing screens representing Figure 24 
used to carry out the simulation in the 3D PE model. The operating frequency is 900 
MHz and same as the previous case. 



Screen Placement (m) 


Figure 27. Side View of two Sinusoidal Hills plus Buildings Profile. 

Figure 28 shows the composite results of the 3D PE model, free-space and Piazzi. 
The solid line is the result of the 3D PE model with 100 m screen width. The dashed line 
is the result of the 50 m screen width. The dashed and dot line is the result of the 25 m 
screen width. The triangle with solid line is the free-space result. And the circles are 
Piazzi’s results. Up to 1500 m away from the base station, the field strength behaves like 
free-space, which corresponds to the peak of the first hill. Then at the first hilltop the 
field behaves like a single knife-edge problem, same as the previous case, where the 
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magnitude of the field strength is 6 dB below the free-space. As in the one hill problem, 
the field does not behave like free-space field when it is close to the base station. The 
field strength minima occur at x= 1000 and 4000 m which do not correspond to the bases 
of the hills as it did in a single round hill case. When we decrease the widths of the 
screen, we see more oscillations as demonstrated in the one hill case. The oscillation 
appears to be more severe in the 25 m screen width case. These oscillations are caused 
by the lateral waves. 



Figure 28. Two Sinusoidal Hills with Rows of Buildings of Uniform Heights and 

Spacing Field Strength. 

We have thus demonstrated that the 3D PE model can support both flat and hilly 
terrains with multiple absorbing screens of uniform heights, equal spacing, and variable 
screen widths. 
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V. CONCLUSIONS AND RECOMMENDATIONS 


A. CONCLUSIONS 

A 3D model based on the parabolic equation (PE) for predicting propagation path 
loss has been developed in this thesis. Two types of terrains were considered: the flat 
earth and the hilly terrain. For the flat earth case, four test problems were examined. The 
results from these cases were compared with results available from the two-ray model, 
the four-ray model, the uniform theory of diffraction (UTD), Lee’s formulation, and the 
numerical integration technique as proposed by Bertoni. For the hilly terrain case, two 
test problems were considered. The results from this case were compared with the results 
presented by Piazzi and Bertoni. For all these test problems, ground was assumed to be 
non-perfectly reflecting, and buildings as perfectly absorbing because all available results 
that were used in the comparisons made these assumptions. 

For the flat earth case, the first test problem was to simulate and determine the 
propagation factor F at the final range over flat earth without any obstacles between the 
transmitter and the receiver. The transmitting antenna had a height FI, = 30 m, and the 
operating frequency was 1 GHz. The receiver located 1000 m from the transmitter. The 
results of the 3D PE model compared excellently with those of the two-ray model. 

The second test problem was to simulate and determine F with a single absorbing 
screen placed between the transmitter and the receiver that had a height Fl^ = 50 m, and a 
width Wk= 49.5 m. The screen representing a building was located 125 m from the 
transmitter and 375 m to the receiver. The transmitting antenna height FI, = 60 m, and the 
operating frequency was 1 GHz. The results of the 3D PE model in this case were 
compared with the results of the four-ray model. Once again, the results had excellent 
agreement. F approached 0 dB as the receiving antenna height increased. In this case, 
both the 3D PE model and the four-ray results approached 0 dB at the receiving antenna 
height of 80 m. 

In the third test problem, nine absorbing screens of uniform heights and equal 
spacing were placed between the transmitter and the receiver. This test case was 
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equivalent to a radiowave propagating over a row of houses in an urban area. The 
screens had a height = 10 m and the widths, Wk- The three different screen widths Wk 
were 12.5 m, 25 m, 50 m. Each screen had a height of 10 m and a spacing distance of 50 
m. Their effects on the propagation factor F were studied. The propagation factor in this 
case was determined at the rooftops. The results of the 3D PE model were compared 
with the UTD results presented by Andersen and the theoretical model results proposed 
by Eee. The 3D PE model results with the screen widths of 25 m and 50 m had excellent 
agreement with the UTD, and both the UTD and the 3D PE model results followed the 
theoretical results only up to the seventh screen. Then up to the ninth screen they were 
0.5 dB higher than the theoretical predictions. Eor the 12.5 m case, the results were 
slightly higher than all cases up to the seventh screen; then they were lower than all cases 
up to the ninth screen. These differences were caused by the lateral waves due to the 
finite screen width. 

The last test problem for the flat earth case involved 120 absorbing screens of 
uniform heights and equal spacing between the transmitter and the receiver. These 
screens represented row of buildings or houses. The transmitting antenna had a height of 
Fit = 125 m, and the operating frequency was 900 MHz. The screens had the height, Hk = 
20 m, and the width, Wk = 50 m. They were spaced 50 m apart. The propagation factor 
was determined at the final range. The results were compared with the results of the 
numerical integration technique presented by Bertoni and the UTD. In the shadow 
region, the results of the 3D PE and the numerical integration method agreed, but the 
result of the UTD was slightly higher, which implied that the UTD overestimated the 
field amplitude. Above the shadow region, the numerical integration method and the 
UTD had excellent agreement, and although the results of the 3D PE model were 0.25 dB 
lower than those two models, nonetheless there was reasonably good agreement with 
them. 

Eor the hilly terrain case, the first test problem was a single rounded hill with 

multiple absorbing screens of uniform heights and equal spacing. The screens had the 

heights Hk = I m. The field strengths were determined on rooftops, and were compared 

with the numerical integration technique results presented by Piazzi and Bertoni. The 

second test problem, the propagation took placed over two hills of sinusoidal shape with 
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multiple absorbing screens of uniform heights and equal spacing. Again, the results were 
compared with the results presented by Piazzi and Bertoni. The screens in this case also 
had the heights Hk = l m In both cases, three variable screen widths Wk were considered: 
25 m, 50 m and 100 m. Their effects on the field strengths were studied. For both cases, 
the transmitting antenna height was Ht = 57 m, and the operating frequency was 900 
MHz. 

For a single round hill, the result of the 3D PE model with screen width of 100 m 
compared excellently with the Piazzi’s result. Before the hilltop, the field strength 
behaved like free-space. Then at the hilltop the field behaved as in a single knife-edge 
case, where the magnitude of the field strength was 6 dB below the free-space. On the 
backside of the hill, the field strength decayed rapidly and linearly with distance, but 
began to increase rapidly to maximum and then decreased slowly with distance. For the 
cases of 25 m and 50 m screen width, the field strengths agreed with thelOO m case and 
Piazzi’s result until after the base of the hill on the backside after which it began to 
increase rapidly with oscillations. These oscillations were caused by the laterally 
diffracted waves. The oscillations were more prominent in the 25 m case than the 50 m 
case. 

For two sinusoidal hills, the field strengths behaved like free-space up to the 
first hilltop similar to the previous case. Then at the first hilltop the fields behaved like a 
single knife-edge problem, where the magnitude of the field strength is 6 dB below the 
free-space. The field strength minima occurred at x = 1000 m and 4000 m which did not 
correspond to the bases of the hills. When the widths of the screen were decreased, the 
oscillations showed up similarly to the single rounded hill case. The oscillations 
appeared to be more severe in the screen width of 25 m. 

Finally, the 3D PE model for predicting propagation path loss in urban areas on 
flat and hilly terrains was developed. Six different test problems were considered. The 
results were compared with the results available in open literature and showed excellent 
agreement. We have also demonstrated that the 3D PE model can support both flat and 
hilly terrains with multiple absorbing screens of uniform heights, equal spacing and 
variable screen widths. 
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B. RECOMMENDATIONS 


Since realistic environments will have non-perfectly reflecting obstacles, a more 
complete 3D PE formulation for non-perfectly reflecting obstacles should be studied as 
well as the effects of the antenna depolarization. In addition, it is recommended that the 
algorithm be implemented in C or Fortran for faster computation, along with the 
automating terrain and building inputs for the model. It may be desirable to add the 
capability to utilize the Digital Terrain Elevation Data (DTED) and Digital Feature 
Attribute Data (DEAD) available at the National Imagery and Mapping Agency (NIMA). 
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APPENDIX A: 3D PE ALGORITHM AND 2RAY MODEL 


Below is the MatLab program listing of the 3D parabolic equation model and the 
2Ray model. This program will compute the propagation factors for the 3D PE with the 
2Ray model over perfectly flat earth and constant admittance plane. The operating 
frequency is 1 GHz. 


o 

% Declaring Contants 

^kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

o 


hold off; 
tic 
c=3e8; 
f=le9; 

Lamda=c/f; 

ko=(2*pi)/Lamda; 

S=2 5; 

Er=5 ; 


% Speed of light 
% Operating frequency 
% Wavelength 
% Wavenumber 

% Ground conductivity in mS/m 
% Relative dielectric constant 


Dz=Lamda; 

Dy=Lamda; 

sigmaz=Lamda; 


Delta Z 
Delta Y 
Current 


in meter 
in meter 

source standard deviation 


Q-kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

o 

% Input Parameters 

^kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

o 


Ht=30; % Transmitting antenna height 


D=1000; 

o, 

0 

Range from 

the 

transmitter 

Dx=100; 

o, 

0 

Incremental 

range 

(Dx) in meter 

Ny=1024; 

o, 

0 

Sample size 

in 

the 

y-direction 

Nz=1024; 

o, 

0 

Sample size 

in 

the 

z-direction 


Q-kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

o 

% Define range in y-direction and z-direction 

^kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

o 


ymax=Ny*Dy/2; 
yl=0:Dy:ymax; 
y=[yl zeros(l,Ny/2-l)]; 


% maximum range in y-direction 

% First half of y range incrementally increase 
% y array 


zmax = Nz*Dz/2; 

zl=0:Dz:zmax; 

z=[zl zeros(l,Nz/2-l)]; 


% Find maximum z range (Zmax) 

% First half of the range of Zmax+1 
% Construct a full z array 
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Q'ifi^iri^iririfirifirifiriririfirifiriririri^ifi!:iri!:iri^ifi^iri!:ifi!:iri^iri^ir-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 

o 

% Define Wavenumbers in Spatial Domain 

^ifififififififififififififififififififififif-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 

o 


kymax=pi/Dy; 

Dky=2 *kymax/Ny; 
kyl=-kymax:Dky:kymax; 
ky=[kyl(:,(Ny/2)+1:Ny+1) kyl 


% Maximum wavenumber in y-direction 
% detla ky 
% Range of ky 
,2:Ny/2)]; 


kzmax=pi/Dz; 

Dkz=2 *kzmax/Nz; 
kzl=-kzmax:Dkz:kzmax; 
kz=[kzl(:,(Nz/2)+1:Nz+1) kzl 


% Maximum wavenumber in z-direction 
% detla kz 
% Range of kz 
,2 :Nz/2)]; 


o 

% Compute Wavenumbers in Frequency Domain, kx 

^-'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 


Ky=meshgrid(ky, 1 

:Nz) ; 

o, 

0 

Create a NyxNz matrix 

of 

ky row-repeat 

Kz=meshgrid(kz, 1 
repeat 

:Ny).'; 

o 

0 

Create a NyxNz matrix 

of 

kz column- 

kx=sqrt(ko^2-Ky. 

^2-Kz.^2); 

o, 

0 

Compute theNy x Nz kx 

mat 

rix 

clear Ky; 


o, 

0 

Clear matrices Ky and 

Kz 

from memory 

clear Kz; 


o, 

0 

To make the algorithm 

run 

faster 


^■k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 

o 

% Compute reflection coefficient 

^kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

o 


Erc=Er+i*18*S/(f/le6); 
Zs=l/sqrt(Ere); 

Gamma=(kz-ko*Zs)./(kz+ko*Zs); 
%gamma=-ones(Nz,Ny); 
gamma=meshgrid(Gamma,1:Ny).'; 


% Complex dielectric 
% Impedance 

% Complex reflection coefficient 
% Setting gamma = -1 

% Create reflection coefficient matrix 
% by column repeating 


Q-kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

o 

% Define the Hanning Window 

^kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

o 


Hya=[]; % 

for t=0:Ny/2; % 

if (t>=0 & t<=3*Ny/8) % 

hy=l; 

elseif (t>=3*Ny/8 & t<=Ny/2)% 
hy=(sin(4*pi*t/Ny))^2; 

end 

Hya=[Hya hy] ; % 

end % 

Yy=fliplr(Hya(:,2:Ny/2 ) ) ; % 

second element to 

o, 

o 

HY=[Hya Yy] ; % 


An empty vector 
Define Ny/2 + 1 points 
Define first 3Ny/8 +1 points 

Define next Ny/8 points 

Construct the first half of Hanning 
window of l+Ny/2 elements 

Flip left to right of the Hanning 
% windows from the 

Ny/2 element for total of [lx(Ny/2-l)] 
Hanning window in y-direction 
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Hzb= [ ]; % 

for t=0:Nz/2; % 

if (t>=0 & t<=3*Nz/8) % 

h z = 1 ; 

elseif (t>=3*Nz/8 & t<=Nz/2)% 
hz=(sin(4*pi*t/Nz))^2; 

end 

Hzb=[Hzb hz]; % 

end % 


An empty vector 
Define Nz/2 + 1 points 
Define first 3Nz/8 tl points 

Define next Nz/8 points 

Construct the first half of Hanning 
window of l+Nz/2 elements 


Yz = fliplr(Hzb(:,2:Nz/2) ) ; 


HZ=[Hzb Yz]'; 


% Flip left to right of the Hanning 
% windows from the second element to 
% Nz/2 element for total of [lx(Nz/2-l)] 
% Hanning window in z-direction 


Hmy=meshgrid(HY,1:Nz) ; 
Hmz=meshgrid(HZ,1:Ny) 


% Row repeat 
% Column repeat 


H3D=(Hmy.*Hmz); 


% 3D Hanning window 


o 

% Gaussian Current Source 

o 


f = (1/ (sqrt (2*pi) *sigmaz) ) *exp (- (z-Ht) .'' 2 / (2*sigmaz^2) ) ; 
f_tilda=Dz*fft(f); % Fourier transform of f(z) 

^■k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 

o 

% Initial H-field at x=0 

^kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

o 

g_tilda=exp (-kz .''2*sigmaz^2/2) ; % initial g tilda 

hyeO_tilda=g_tilda.*cos(kz*Ht); % Initial hye_tilda(0+,ky,kz) field 

hyoO_tilda=-i*g_tilda.*sin(kz*Ht); % Initial hyo_tilda(0+,ky,kz) field 

% Initial even H-field, column repeat 
HyeO_tilda=meshgrid(hyeO_tilda,1:Ny).'; 

% Initial odd H-field, column repeat 
HyoO_tilda=meshgrid(hyoO_tilda,1:Ny).'; 

% Include the reflection coefficient 
Hy0_tilda=0.5*(1+gamma).*Hye0_tilda+0.5*(1-gamma).*HyoO_tilda; 

Hy_tilda=HyO_tilda.*H3D; % Hy_tilda(0+,ky,kz) at x=0 

Exp2=exp(i*kx*Dx); % Marching range 
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o 


% 3D Parabolic Equation Basic Algorithm 

^ififififififififififififififififififififififif-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 


for x=0:Dx:D-Dx; 

Hy_tilda_Dx=Hy_tilda.*Exp2; 


% Hy_tilda(Dx, ky, kz ) 


Hy=Dkz*Dky*(Ny*Nz*ifft2(Hy_tilda_Dx))/(2*pi)^2; 

% ifft2 wrt to ky and kz 


Hy_H=Hy.*H3D; 

Hyz1=Hy_H(l:Nz/2 + l, :) ; 

HyzO = flipud(Hy_H(2:Nz/2, :) ) ; 

Hye=[Hyzl; HyzO]; 

Hyo=[Hyzl; -HyzO]; 

Hye_tilda=(Dz*Dy)*fft2(Hye) ; 
Hye_tilda_H=Hye_tilda.*H3D; 

Hyo_tilda=(Dz*Dy)*fft2(Hyo) ; 
Hyo_tilda_H=Hyo_tilda.*H3D; 


% Apply the Hamming window 
% H-field for z > 0 
% H-field for z < 0 
% Even H-field 
% Odd H-field 

% Take the EFT of the even H~ 

% Apply the Hanning window in freq 

% Take the EFT of the odd H~ 

% Apply the Hanning window in freq 


% Multiply by the reflection coefficient 
Hye_tilda_g=0.5*(1+gamma).*Hye_tilda_H; 

% Multiply by the reflection coefficient 
Hyo_tilda_g=0.5*(1-gamma).*Hyo_tilda_H; 

Hy_tildal=Hye_tilda_g + Hyo_tilda_g; % Hy_tilda(x,ky,kz) 
Hy_tilda=Hy_tildal.*H3D; % Apply the Hanning window 


end 

o 

% Compute Propagation Factor at Final Range 

^■k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 

o 

J=2; % Pick a column of H-field at final range 

rho=sqrt(D^2+(J)^2); % Distance from transmitter to receiver 

R=sqrt( (z ( :,1:Nz/2 + 1)-Ht) .^2 + rho^2); 

% Distance from transmitter to receiver 
theta=asin(rho./R); % Angle Theta in spherical coordinate 

phi=asin(J/rho) ; % Angle Phi in spherical coordinate 

g_tilda_fs=exp(-((cos(theta)*ko).^2*sigmaz^2)/2); 

% Free-space g_tilda 

Hyfs=-i*ko*sin(theta)*cos(phi).*g_tilda_fs.*exp(i*ko*R)./(4*pi*R); 

% Free-space H-field 

A=Hy_H(1:Nz/2+1,2); % z > 0 of a column of the H-field 

F=20*logl0(abs(A./Hyfs')); % Propagation factor 

save H:\THESIS_2001\MATLAB_CODES\PERF2Ray F -ASCII 
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Q'ifi^iri^iririfirifirifiriririfiriririfiriririfi^iri^iri^ifi^iri!:ifi^iri^iri^iri!:-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 

o 


% 2Ray-Model 

^ifififififififififififififififififififififif-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 

o 


J=2; 

o, 

0 

rho=sqrt(D^2+(J)^2); 

o, 

0 

receiver 


Hr=z(1,1:Nz/2+l); 

o, 

0 

DR=2*Ht*Hr/rho; 

o 

0 

o. 

psi=atan((Ht+Hr)/rho); 

0 

o, 

0 


A column final H-field 

Distance from transmitter to the 

The receiving antenna height 
Path difference between direct and 
ground reflected waves 
Grazing angle 


Z2ray=sqrt(Erc-cos(psi).^2)/Ere; % Normalized surface impedance Z(psi) 
gamma2ray=((sin(psi)-Z2ray)./(sin(psi)+Z2ray)); % Reflection coefficient 
Tworay=abs(l+gamma2ray.*exp(i*ko*DR)); % 2ray model propagation factor 
F2ray=20*logl0(Tworay) ; % 2ray model propagation factor in dB 

save H:\THESIS_2001\Base_Line_Codes\RF2Ray F2ray -ASCII 
save H:\THESIS_2001\Base_Line_Codes\3DHanning H3D -ASCII 

9''k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

* Plot results 

^-'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 


figure(1) 
mesh(H3D) 
view(-37.5, 70) 

figure (2) 

y=0:Dz:(Nz)*Dz/2; 

plot(F,y,F2ray,y, 'o-r') 

axis ( [-30 6 0 Nz*Dz/2]) 

legend( '3D Parabolic' , '2Ray' ) 

toe 
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APPENDIX B: A 3D PE MODEL AND 4RAY MODEL 


Below is the MatLab program listing of the 3D parabolic equation model and the 
4Ray model. This program will compute the propagation factors for the 3D PE model 
and the 4Ray model over perfectly flat earth and constant admittance plane. A single 
absorbing knife-edge is placed between he transmitter and receiver. The operating 
frequency is 1 GHz. 


o 

% Declaring Contants 

o 


hold off; 
tic 
c=3e8; 
f=le9; 

Lamda=c/f; 

ko=(2*pi)/Lamda; 

S=400; 

Er=10; 


% Speed of light 
% Operating frequency 
% Wavelength 
% Wavenumber 

% Ground conductivity in mS/m 
% Relative dielectric constant 


Dz=Lamda; 

Dy=Lamda; 

sigmaz=Lamda; 


% Delta Z in meter 
% Delta Y in meter 
% Standard deviation of Z 


o 

% Input Parameters 

^■k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 

o 


Ht=60; 
Hk=5 0 ; 
Wk=49.8; 


% Transmitting antenna height 
% Knife edge height 
% Knife edge width 


dl=125; 
d2=375; 
D=500; 
Dx=125; 


% Distance from transmitter to knife 
% Distance from knife edge to receiver 
% Range from the transmitter 
% Incremental range (Dx) in meter 


Ny=1024; 

Nz=1024; 


% Sample size in the y-direction 
% Sample size in the z-direction 


Q-kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

o 

% Define range in y-direction and z-direction 

^kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

o 


ymax=Ny*Dy/2; % 
yl=0:Dy:ymax; % 
y=[yl zeros(1,Ny/2-1)]; % 


maximum range in y-direction 
First half of y range 
y array 


zmax = Nz*Dz/2; 

zl=0:Dz:zmax; 

z=[zl zeros(l,Nz/2-l)]; 


% Find maximum z range (Zmax) 

% First half of the range of Zmax+1 
% Construct a full z array 
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o 

% Define Wavenumbers in Spatial Domain 

^kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

o 


kymax=pi/Dy; 

Dky=2 *kymax/Ny; 
kyl=-kymax:Dky:kymax; 
ky=[kyl(:,(Ny/2)+1:Ny+1) kyl 


% Maximum wavenumber in y-direction 
% detla ky 
% Range of ky 
,2:Ny/2)]; 


kzmax=pi/Dz; 

Dkz=2 *kzmax/Nz; 
kzl=-kzmax:Dkz:kzmax; 
kz=[kzl(:,(Nz/2)+1:Nz+1) kzl 


% Maximum wavenumber in z-direction 
% detla kz 
% Range of kz 
,2:Nz/2)]; 


Q-kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

o 

% Compute Wavenumbers in Frequency Domain, kx 

^kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

o 


Ky=meshgrid(ky,1:Nz); 
Kz=meshgrid(kz,1:Ny) . ' ; 
repeat 

kx=sqrt(ko^2-Ky.^2-Kz.^2) ; 


% Create a NyxNz matrix 
% Create a NyxNz matrix 


of ky row-repeat 
of kz column- 


% Compute theNy x Nz kx matrix 


clear Ky; 
clear Kz; 


% Clear matrices Ky and Kz from memory 
% To make the algorithm run faster 


o 

% Compute reflection coefficient 

^-kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

o 


Erc=Er+i*18*S/(f/le6); 
Zs=l/sqrt(Ere); 

Gamma=(kz-ko*Zs)./(kz+ko*Zs); 
%gamma=-ones(Nz,Ny); 
gamma=meshgrid(Gamma,1:Ny).'; 


% Complex dielectric 
% Impedance 

% Complex reflection coefficient 
% Setting gamma = -1 

% Create reflection coefficient matrix 
% by column repeating 


Q-kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

o 

% Define the Hanning Window 

^kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

o 


Hya=[]; % 

for t=0:Ny/2; % 

if (t>=0 & t<=3*Ny/8) % 

hy=l; 

elseif (t>=3*Ny/8 & t<=Ny/2)% 
hy=(sin(4*pi*t/Ny))^2; 

end 

Hya=[Hya hy] ; % 

end % 

Yy=fliplr(Hya(:,2:Ny/2 ) ) ; % 


An empty vector 
Define Ny/2 + 1 points 
Define first 3Ny/8 +1 points 

Define next Ny/8 points 

Construct the first half of Hanning 
window of l+Ny/2 elements 

Flip left 


second element to 
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to right of the Hanning 
% windows from the 


HY=[Hya Yy]; 


% Ny/2 element for total of [lx(Ny/2-l)] 
% Hanning window in y-direction 


Hzb= [ ]; % 

for t=0:Nz/2; % 

if (t>=0 & t<=3*Nz/8) % 

h z = 1 ; 

elseif (t>=3*Nz/8 & t<=Nz/2)% 
hz=(sin(4*pi*t/Nz) ) '' 2 ; 

end 

Hzb=[Hzb hz]; % 

end % 


An empty vector 
Define Nz/2 + 1 points 
Define first 3Nz/8 tl points 

Define next Nz/8 points 


Construct the first half of Hanning 
window of l+Nz/2 elements 


Yz = fliplr(Hzb(:,2:Nz/2) ) ; 


HZ=[Hzb Yz]'; 


% Flip left to right of the Hanning 
% windows from the second element to 
% Nz/2 element for total of [lx(Nz/2-l)] 
% Hanning window in z-direction 


Hmy=meshgrid(HY,1:Nz) ; 
Hmz=meshgrid(HZ,1;Ny) 


% Row repeat 
% Column repeat 


H3D=(Hmy.*Hmz); 


% 3D Hanning window 


o 

% Gaussian Current Source 

o 


f = (1/ (sqrt (2*pi) *sigmaz) ) *exp (- (z-Ht) . "'2/ (2*sigmaz''2) ) ; 
f_tilda=Dz*fft(f); % Fourier transform of f(z) 

^■k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 

o 

% Initial H-field at x=0 

^kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

o 


g_tilda=exp (-kz .''2*sigmaz^2/2) ; % initial g tilda 

hyeO_tilda=g_tilda.*cos(kz*Ht); % Initial hye_tilda(0+,ky,kz) field 

hyoO_tilda=-i*g_tilda.*sin(kz*Ht); % Initial hyo_tilda(0+,ky,kz) field 

% Initial even H-field, column repeat 
HyeO_tilda=meshgrid(hyeO_tilda,1:Ny).'; 

% Initial odd H-field, column repeat 
HyoO_tilda=meshgrid(hyoO_tilda,1:Ny).'; 

% Include the reflection coefficient 
Hy0_tilda=0.5*(1+gamma).*Hye0_tilda+0.5*(1-gamma).*HyoO_tilda; 

Hy_tilda=HyO_tilda.*H3D; % Hy_tilda(Ot,ky,kz) at x=0 

Exp2=exp(i*kx*Dx); % Marching range 
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o 

% 3D Parabolic Equation Basic Algorithm 

^ifififififififififififififififififififififif-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 

o 

C=0; % Setting the counter 


for x=0:Dx:D-Dx; 

Hy_tilda_Dx=Hy_tilda.*Exp2; 

% Hy_tilda(Dx, ky, kz ) 


Hy=Dkz*Dky*(Ny*Nz*ifft2(Hy_tilda_Dx))/(2*pi)^2; % ifft2 wrt to ky and 


kz 


Hy_H=Hy.*H3D; % Apply 

the Hamming window in spatial domain 

C=C+Dx; 

% Counter 


if (C == dl) % When X= 125 M zero out the field 

Hy_H(1:round(Hk/Dz),1:round(Wk/Dy/2))=0; 

Hy_H(1:round(Hk/Dz),(Ny-round(Wk/Dy/2));Ny)=0; 

end 


Hyz1=Hy_H(l:Nz/2 + l, :) ; 

% H-field for z > 0 

HyzO = flipud(Hy_H(2:Nz/2, :) ) ; 

% H-field for z < 0 

Hye=[Hyzl; HyzO]; 

% Even H-field 

Hyo=[Hyzl; -HyzO]; 

% Odd H-field 

Hye_tilda=(Dz*Dy)*fft2(Hye) ; 
Hye_tilda_H=Hye_tilda.*H3D; 

% Take the FFT of the even H~ 

% Apply the Hanning window in freq 

Hyo_tilda=(Dz*Dy)*fft2(Hyo) ; 
Hyo_tilda_H=Hyo_tilda.*H3D; 

o, 

0 

% Take the FFT of the odd H~ 

% Apply the Hanning window in freq 
Multiply by the reflection coefficient 


Hye_tilda_g=0.5*(1+gamma).*Hye_tilda_H; 

% Multiply by the reflection coefficient 
Hyo_tilda_g=0.5*(1-gamma).*Hyo_tilda_H; 

Hy_tildal=Hye_tilda_g + Hyo_tilda_g; % Hy_tilda(x,ky,kz) 


Hy_tilda=Hy_tildal.*H3D; 

% Apply the Hanning window 

end 



o 

% Compute the propagation factor at the final range 

^-'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

J=2; % First column final H-field 

rho=sqrt(D^2+(J)^2); % Distance from Tx base of Rx 

R=sqrt ( ( z ( : , 1: Nz/2 + 1)-Ht) . "'2 + rho''2); % Distance from Tx to Rx 


theta=asin(rho./R) ; 
phi=asin(J/rho) ; 

g_tilda_fs=exp(-((cos(theta)*ko) 

% Angle Theta 
% Angle Phi 

1 .^2*sigmaz^2)/2); % Free-space g_tilda 


Hyfs=-i*ko*sin(theta)*cos(phi).*g_tilda_fs.*exp(i*ko*R)./(4*pi*R); % 
Free-space H-field 
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A=Hy_H(1:Nz/2 + l, 2) ; 
F=20*logl0(abs(A./Hyfs')) ; 


1st column of the final H-field 
Propagation factor 


o, 

o 

o, 

o 


save H:\THESIS_2001\MATLAB_CODES\PERF4Ray F -ASCII 

o 

% Four-Ray Model 

o 


% The radius of any member n of the 
perpendicular to the LOS path 

Hl=sqrt((Lamda*dl*d2)/(dl+d2)); 

hl=Hk-(Ht*d2+z(:,1:Nz/2+l)*dl)/D; 
Gamma1=1; 
xl=sqrt(2)*hl/Hl; 

Cl=mfun( 'FresnelC' ,xl); 

Sl=mfun( 'Fresnels' ,xl); 

Fnl=Cl-j*Sl; 

Fl=0.5*(l-(l+j)*Fnl)*Gammal; 
rl=sqrt(D^2+(Ht-z(:,l:Nz/2 + l)) .^2) ; 
E01=exp(-j*ko*rl)./rl; 

Edl=E01.*F1; 

h2=Hk+(Ht*d2-z(:,1:Nz/2+l)*dl)/D; 
psiA=atan((Hk+Ht)/dl); 

%GammaA=(sin(psiA)-Zs)/(sin(psiA)+Z 
GammaA=-l; 
x2=sqrt(2)*h2/Hl; 

C2=mfun( 'FresnelC' ,x2); 

S2=mfun( 'FresnelS' ,x2) ; 

Fn2=C2-j*S2; 

F2=0.5*(l-(l+j)*Fn2)*GammaA; 

r2 = sqrt(D^2+(Ht + z(:,l:Nz/2 + l) ) .^2) ; 

E02=exp(-j*ko*r2)./r2; 

Ed2=E02.*F2; 

h3=Hk-(Ht*d2-z(:,1:Nz/2+l)*dl)/D; 
psiB=atan((Hk+z(:,l:Nz/2+l))/d2); 

% GammaB=(sin(psiB)-Zs)./(sin(psiB) 
GammaB=-l; 
x3=sqrt(2)*h3/Hl; 

C3=mfun( 'FresnelC' ,x3); 

S3=mfun ('FresnelS',x3); 

Fn3=C3-j*S3; 

F3=0.5*(l-(l+j)*Fn3).*GammaB; 
r3=sqrt(D^2+(Ht+z(:,l:Nz/2+l)).^2); 
E03=exp(-j*ko*r3)./r3; 

Ed3=E03.*F3; 

h4=Hk+(Ht*d2+z(:,1:Nz/2+l)*dl)/D; 

GammaC=GammaA*GammaB; 

x4=sqrt(2)*h4/Hl; 

C4=mfun( 'FresnelC' ,x4); 

S4=mfun( 'FresnelS' ,x4) ; 

Fn4=C4-j*S4; 


family of ellipsoids in a plane 


% 1st clearance height 
% Reflection coefficient for case 1 

% Real part of the Fresnel Integral 
% Imag part of the Fresnel Integral 
% Solution for the Fresnel Integral 
% First knife-edge normalized field 


% 2nd clearance height 
% Grazing angle of case 2 
s);% Refl. coefficient for case 2 


% Real part of the Fresnel Integral 
% Imag part of the Fresnel Integral 
% Solution for the Fresnel Integral 
% 2nd knife-edge normalized field 


% 3rd clearance height 
% Grazing angle of case 3 
+Zs); % Refl. coefficient for case 3 

% Real part of the Fresnel Integral 
% Imag part of the Fresnel Integral 
% Solution for the Fresnel Integral 
% 3rd knife-edge normalized field 


% 4th clearance height 
% Reflection coefficient for case 4 

% Real part of the Fresnel Integral 
% Imag part of the Fresnel Integral 
% Solution for the Fresnel Integral 
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F4=0.5*(l-(l+j)*Fn4).*GammaC; 
r4=sqrt(D^2+(Ht-z(:,l:Nz/2+l)).^2); 
E04=exp(-j*ko*r4)./r4; 

Ed4=E04.*F4; 


4th knife-edge normalized field 


Fknife=(Edl+Ed2+Ed3+Ed4)./EOl; 


Total normalized fields 


FknifedB=20*logl0(abs(Fknife)); 


Total normalized fields in dB 


save H:\THESIS_2001\MATLAB_CODES\RF4Ray FknifedB -ASCII 



Plot results 



figure(1) 

y=0:Dz:(Nz)*Dz/2; 

plot(F,y,FknifedB,y, 'o-r' ) 

axis ( [-30 6 0 Nz*Dz/2]) 

%xlabel('Relative Propagation Factor F(dB)'),ylabel('Altitude(m)') 
%title(['Ny = 'num2str(Ny) ' Nz = ' num2str(Nz) ' Dx = 'num2str(Dx) ' 
meters' ' Ht = ' num2str(Ht) ' meters']) 
legend( 'Parabolic' , 'KnifedB' ) 
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APPENDIX C: A 3D PE MODEL FOR THE MULTIPLE SCREENS 
OF UNIFORM HEIGHTS AND EQUAL SPACING 


Below is the MatLab program listing of the 3D parabolic equation model use to 
compute the propagation factors with multiple absorbing screens of uniform heights and 
equal spacing over perfectly flat earth and constant admittance plane. The operating 


frequency is 1 GHz and 900 MHz. 


^ifififif-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 

o 

% Declaring Contants 


hold off; 

tic 

c=3e8; 

f=0.9e9; 

Lamda=c/f; 

ko=(2*pi)/Lamda; 

S=400; 

Er=10; 

Dz=Lamda; 
DY=Lamda; 
sigmaz=Lamda; 


% Speed of light 
% Operating frequency 
% Wavelength 
% Wavenumber 

% Ground conductivity in mS/m 
% Relative dielectric constant 

% Delta Z in meter 
% Delta Y in meter 

% Current source standard deviation 


^ifi!:ifif-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 

o 

% Input Parameters 


Ht=125; 

Hk=20; 

Wk=5 0; 

NS=120; 

Dns=50 ; 

M=50:Dns:Dns*NS; 

D=Dns*(NS+1); 
Dx=2 5; 


% Transmitting antenna height 

% Knife edge height 
% Knife edge width 

% Number of screens 
% Screen separation 
% Location of the screens 

% Range from the transmitter 
% Incremental range (Dx) in meter 


Ny=1024; 

Nz=1024; 


% Sample size in the y-direction 
% Sample size in the z-direction 


^ifififif-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 

o 

% Define range in y-direction and z-direction 

^■k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 


ymax=Ny*Dy/2; % maximum range in y-direction 

yl=0:Dy:ymax; % First half of y range 

y=[Yl zeros(1,Ny/2-1)]; % y array 
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zmax = Nz*Dz/2; % Find maximum z range (Zmax) 

zl=0:Dz:zmax; % First half of the range of Zmax+1 

z=[zl zeros(1,Nz/2-1)]; % Construct a full z array 


o 

% Define Wavenumbers in Spatial Domain 

^ifififififififififififififififififififififif-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 

o 

kYmax=pi/Dy; % Maximum wavenumber in y-direction 

Dky=2*kymax/Ny; % detla ky 

kyl=-kymax:Dky:kymax; % Range of ky 

ky=[kyl(: , (Ny/2)+1:Ny+l) kyl(:,2:Ny/2)]; 

kzmax=pi/Dz; % Maximum wavenumber in z-direction 

Dkz=2*kzmax/Nz; % detla kz 

kzl=-kzmax:Dkz:kzmax; % Range of kz 

kz=[kzl(:, (Nz/2)+1:Nz + l) kzl(;,2:Nz/2) ] ; 

o 

% Compute Wavenumbers in Frequency Domain, kx 


Ky=meshgrid(ky, 1 

:Nz) ; 

o, 

0 

Create a NyxNz matrix 

of 

ky row-repeat 

Kz=meshgrid(kz, 1 
repeat 

:Ny).'; 

o, 

0 

Create a NyxNz matrix 

of 

kz column- 

kx=sqrt (ko''2-Ky. 

''2-Kz.''2) ; 

o, 

0 

Compute theNy x Nz kx 

mat 

r ix 

clear Ky; 


o, 

0 

Clear matrices Ky and 

Kz 

from memory 

clear Kz; 


o, 

0 

To make the algorithm 

run 

faster 


^ifififif-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 

o 

% Compute reflection coefficient 

^■k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 

o 

Erc=Er+i*18*S/(f/le6); % Complex dielectric 

Zs=l/sqrt(Ere); % Impedance 

Gamma=(kz-ko*Zs)./(kz+ko*Zs); % Complex reflection coefficient 
%gamma=-ones(Nz,Ny) ; % Setting gamma = -1 

gamma=meshgrid(Gamma,1:Ny); % Create reflection coefficient matrix 

% by column repeating 

^kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

o 

% Define the Hanning Window 

^kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 


An empty vector 
Define Ny/2 + 1 points 
Define first 3Ny/8 +1 points 


Hya= [ ] ; 
for t=0:Ny/2; 

if (t>=0 & t<=3*Ny/8) 
hy=l; 

elseif (t>=3*Ny/8 & t<=Ny/2)% Define next Ny/8 points 
hy=(sin(4*pi*t/Ny))^2; 

end 

Hya=[Hya hy]; 

end 


% Construct the first half of Hanning 
% window of l+Ny/2 elements 
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Yy=fliplr(Hya(:,2:Ny/2 ) ) ; 
second element to 
HY=[Hya Yy]; 


% Flip left to right of the Hanning 
% windows from the 

% Ny/2 element for total of [lx(Ny/2-l)] 
% Hanning window in y-direction 


Hzb= [ ]; % 

for t=0:Nz/2; % 

if (t>=0 & t<=3*Nz/8) % 

h z = 1 ; 

elseif (t>=3*Nz/8 & t<=Nz/2)% 
hz=(sin(4*pi*t/Nz))^2; 

end 

Hzb=[Hzb hz]; % 

end % 


An empty vector 
Define Nz/2 + 1 points 
Define first 3Nz/8 tl points 

Define next Nz/8 points 

Construct the first half of Hanning 
window of l+Nz/2 elements 


Yz=fliplr(Hzb(:,2:Nz/2)) ; 


HZ=[Hzb Yz]'; 


% Flip left to right of the Hanning 
% windows from the second element to 
% Nz/2 element for total of [lx(Nz/2-l)] 
% Hanning window in z-direction 


Hmy=meshgrid(HY,1:Nz) ; 
Hmz=meshgrid(HZ,1:Ny) 


% Row repeat 
% Column repeat 


H3D=(Hmy.*Hmz); 


% 3D Hanning window 


o 

% Gaussian Current Source 

o 


f = (1/ (sqrt (2*pi) *sigmaz) ) *exp (- (z-Ht) .'' 2 / (2*sigmaz''2) ) ; 
f_tilda=Dz*fft(f); % Fourier transform of f(z) 

^■k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 

o 

% Initial H-field at x=0 

S'kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

o 


g_tilda=exp (-kz .''2*sigmaz^2/2) ; % initial g tilda 

hyeO_tilda=g_tilda.*cos(kz*Ht); % Initial hye_tilda(0+,ky,kz) field 

hyoO_tilda=-i*g_tilda.*sin(kz*Ht); % Initial hyo_tilda(Ot,ky,kz) field 

% Initial even H-field, column repeat 
HyeO_tilda=meshgrid(hyeO_tilda,1:Ny).'; 

% Initial odd H-field, column repeat 
HyoO_tilda=meshgrid(hyoO_tilda,1:Ny).'; 

% Include the reflection coefficient 
Hy0_tilda=0.5*(1+gamma).*Hye0_tilda+0.5*(1-gamma).*HyoO_tilda; 

Hy_tilda=HyO_tilda.*H3D; % Hy_tilda(Ot,ky,kz) at x=0 

Exp2=exp(i*kx*Dx); % Marching range 
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o 


% 3D Parabolic Equation Basic Algorithm 

^ifififififififififififififififififififififif-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 

o 

C=0; % Setting the counter 


for x=0:Dx:D-Dx; 

Hy_tilda_Dx=Hy_tilda.*Exp2; 

% Hy_tilda(Dx, ky, kz ) 

% ifft2 wrt to ky and kz 


Hy=Dkz*Dky*(Ny*Nz*ifft2(Hy_tilda_Dx))/(2*pi)^2; 

Hy_H=Hy.*H3D; % Apply the Hamming window in spatial domain 

C=C+Dx; % Counter 

while M(M==C) 

Hy_H(1:round(Hk/Dz),1:round(Wk/Dy/2))=0 ; 

Hy_H(1:round(Hk/Dz),(Ny-round(Wk/Dy/2)):Ny)=0; 
break 

end 


Hyz1=Hy_H(l:Nz/2 + l, :) ; 

% H-field for z > 0 

HyzO = flipud(Hy_H(2:Nz/2, :) ) ; 

% H-field for z < 0 

Hye=[Hyzl; HyzO]; 

% Even H-field 

Hyo=[Hyzl; -HyzO]; 

% Odd H-field 

Hye_tilda=(Dz*Dy)*fft2(Hye) ; 
Hye_tilda_H=Hye_tilda.*H3D; 

% Take the EFT of the even H~ 

% Apply the Hanning window in freq 

Hyo_tilda=(Dz*Dy)*fft2(Hyo) ; 
Hyo_tilda_H=Hyo_tilda.*H3D; 

o, 

0 

% Take the EFT of the odd H~ 

% Apply the Hanning window in freq 
Multiply by the reflection coefficient 


Hye_tilda_g=0.5*(1+gamma).*Hye_tilda_H; 

% Multiply by the reflection coefficient 
Hyo_tilda_g=0.5*(1-gamma).*Hyo_tilda_H; 

Hy_tildal=Hye_tilda_g + Hyo_tilda_g; % Hy_tilda(x,ky,kz) 
Hy_tilda=Hy_tildal.*H3D; % Apply the Hanning window 

end 

^■k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 

o 

% Compute The Propagation Factor 

^kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

o 


II 

% First column final H-field 

rho=sqrt(D^2+(J)^2); 

% Distance from Tx to base Rx 


R=sqrt((z(:,1:Nz/2+1)-Ht).^2 + rho^2); % Distance from Tx to Rx 


theta=asin(rho./R); 
phi=asin(J/rho); 

g_tilda_fs=exp(-((cos(theta)*ko) 

% Angle Theta 
% Angle Phi 

1 .^2*sigmaz^2)/2); % Free-space g_tilda 

% Free-space H-field 



Hyfs=-i*ko*sin(theta)*cos(phi).*g_tilda_fs.*exp(i*ko*R)./(4*pi*R); 
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A=Hy_H(1:Nz/2+l,2); % z > 0 of the first column of the final H-field 
F=(abs(A./Hyfs.')); % Propagation factor relative to the free-space 

save H:\THESIS_2001\MATLAB_CODES\propfact F -ASCII 

S-'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

% Plot results 

o 

figure(1) 

%Z=linspace(0,3*Nz/8*Dz, 513) ; 

Z=(0:Dz:Nz/2*Dz) ' ; 

Y=(Z-Hk)/Lamda; 

save H:\THESIS_2001\MATLAB_CODES\recantenna Y -ASCII 
plot(Y,F) 

axis([-30 60 0 1.8 ]) 

ylabel (' Relative Propagation Factor, F '), xlabel (' Receiving Antenna in 
Wavelegth' ) 

title (['Ny = 'num2str(Ny) ' Nz = ' num2str(Nz) ' Dx = 'num2str(Dx) ' 

meters' ' Ht = ' num2str(Ht) ' meters']) 

grid 

toe % Stop stopwatch 
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APPENDIX D: A 3D PE MODEL FOR A SINGLE ROUND HILL 


Below is the MatLab program listing of the 3D parabolic equation model use to 
compute the field strengths for a single round hill with multiple absorbing screens of 
uniform heights and equal spacing. These screens represent a row of buildings or houses 
in residential areas of a well built-up city. The operating frequency is 900 MHz. 

9''k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

% Declaring Contants 

o 


hold off; 

tic 

c=3e8; 

f=0.9e9; 

Lamda=c/f; 

ko=(2*pi)/Lamda; 

S=400; 

Er=10; 


% Speed of light 
% Operating frequency 
% Wavelength 
% Wavenumber 

% Ground conductivity in mS/m 
% Relative dielectric constant 


Dz=Lamda; 

Dy=Lamda; 

sigmaz=Lamda; 


% Delta Z in meter 
% Delta Y in meter 

% Current source standard deviation 


^ififififif-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 


% Input Parameters 

^■k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 


Ht=57; 
Hbs=7; 


% Transmitting antenna height 
% Receiving antenna height above KE's 


Wk=100; 
Dns=50 ; 


% Knife edge width of 25 meter 
% Screen separation 


Dist=10000; % 
xs=-l000:Dns:1000; % 
ys=sqrt(Ie4^2-xs.^2)-(9.95e3) ;% 
xsl=1000+Dns:Dns:Dist+Dns; % 
ysl=zeros(l,(length(xsl))); % 
xs2=[xs xsl]; 

Hk=[ys ysl]; % 
Hk(l)=0; % 
Hk (length(xs))=0 ; 


Final range 
Horizontal hill width 
Equation for circle 

Equal heigth and equal distance screen 
Equal heigth screen height 

Knife edge height of meter 

When screen height less than zeros 


NS=length(xs2); 
M=Dns:Dns:Dns*(NS) ; 


% Number of screens 
% Location of the screens 


D=Dns*(NS); 
Dx=5 0 ; 


% Range from the transmitter 
% Incremental range (Dx) in meter 


Ny=1024; 

Nz=1024; 


% Sample size in the y-direction 
% Sample size in the z-direction 
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% Define range in y-direction and z-direction 

^ifififififififififififififififififififififif-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 

o 


ymax=Ny*Dy/2; 
yl=0:Dy:ymax; 
y=[yl zeros(1,Ny/2-1) 

zmax = Nz*Dz/2; 

zl=0:Dz:zmax; 

z=[zl zeros(1,Nz/2-1) 


% maximum range 
% First half of 
]; % y array 

% Find maximum 
% First half of 
] ; % Construct a f 


in y-direction 
y range 

range (Zmax) 
the range of Zmax+1 
11 z array 


o 

% Define Wavenumbers in Spatial Domain 

^'’k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 


kymax=pi/Dy; 

Dky=2 *kymax/Ny; 
kyl=-kymax:Dky:kymax; 
ky=[kyl(:,(Ny/2)+1:Ny+1) kyl( 

kzmax=pi/Dz; 

Dkz=2 *kzmax/Nz; 
kzl=-kzmax:Dkz:kzmax; 
kz=[kzl(:, (Nz/2)+1:Nz + 1) kzl ( 


% Maximum wavenumber in y-direction 
% detla ky 
% Range of ky 
,2:Ny/2)]; 

% Maximum wavenumber in z-direction 
% detla kz 
% Range of kz 
,2:Nz/2)]; 


o 

% Compute Wavenumbers in Frequency Domain, kx 

o 


Ky=meshgrid(ky,1:Nz) ; 
Kz=meshgrid(kz,1:Ny) . ' ; 
repeat 

kx=sqrt (ko-'C-Ky. ''2-Kz . ''2) ; 


% Create 
% Create 


a NyxNz matrix 
a NyxNz matrix 


of ky row-repeat 
of kz column- 


% Compute theNy x Nz kx matrix 


clear Ky; 
clear Kz; 


% Clear matrices Ky and Kz from memory 
% To make the algorithm run faster 


Q'iri^iri^ifirifi^iririfiriririfi^ifiriririririfi^ifiriririfiriri!:ifi^iri^iri^if-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 

o 

% Compute reflection coefficient 

^ififif-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 

o 


Erc=Er+i*18*S/(f/le6); 
Zs=l/sqrt(Ere); 

Gamma=(kz-ko*Zs)./(kz+ko*Zs); 
%gamma=-ones(Nz,Ny); 
gamma=meshgrid(Gamma,1:Ny).'; 


% Complex dielectric 
% Impedance 

% Complex reflection coefficient 
% Setting gamma = -1 

% Create reflection coefficient matrix 
% by column repeating 


o 

% Define the Hanning Window 

S'if'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 


Hya= [ ] ; 

for t=0:Ny/2; 


% An empty vector 
% Define Ny/2 + 1 points 
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if (t>=0 & t<=3*Ny/8) 
hy=l; 

elseif (t>=3*Ny/8 & t<=Ny/2 
hy=(sin(4*pi*t/Ny))^2; 

end 

Hya=[Hya hy]; 

end 


% Define first 3Ny/8 +1 points 
% Define next Ny/8 points 

% Construct the first half of Hanning 
% window of l+Ny/2 elements 


Yy=fliplr(Hya(:,2:Ny/2)) ; 
second element to 
HY=[Hya Yy]; 


% Flip left to right of the Hanning 
% windows from the 

% Ny/2 element for total of [lx(Ny/2-l)] 
% Hanning window in y-direction 


Hzb= [ ]; % 

for t=0:Nz/2; % 

if (t>=0 & t<=3*Nz/8) % 

h z = 1 ; 

elseif (t>=3*Nz/8 & t<=Nz/2)% 
hz=(sin(4*pi*t/Nz))^2; 

end 

Hzb=[Hzb hz]; % 

end % 


An empty vector 
Define Nz/2 + 1 points 
Define first 3Nz/8 +1 points 

Define next Nz/8 points 

Construct the first half of Hanning 
window of l+Nz/2 elements 


Yz = fliplr(Hzb(:,2:Nz/2) ) ; 


HZ=[Hzb Yz]'; 


% Flip left to right of the Hanning 
% windows from the second element to 
% Nz/2 element for total of [lx(Nz/2-l)] 
% Hanning window in z-direction 


Hmy=meshgrid(HY,1:Nz) ; 
Hmz=meshgrid(HZ,1:Ny) ' ; 


% Row repeat 
% Column repeat 


H3D=(Hmy.*Hmz); 


% 3D Hanning window 


o 

% Gaussian Current Source 

S-'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 


f = (1/ (sqrt (2*pi) *sigmaz) ) *exp (- (z-Ht) .'' 2 / (2*sigmaz^2) ) ; 
f_tilda=Dz*fft(f); % Fourier transform of f(z) 

^ifififirifirifir-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 

o 

% Initial H-field at x=0 

^-'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 


g_tilda=exp(-kz.^2*sigmaz^2/2); % initial g tilda 

hyeO_tilda=g_tilda.*cos(kz*Ht); % Initial hye_tilda(0+,ky,kz) field 

hyoO_tilda=-i*g_tilda.*sin(kz*Ht); % Initial hyo_tilda(0+,ky,kz) field 

% Initial even H-field, column repeat 
HyeO_tilda=meshgrid(hyeO_tilda,1:Ny).'; 

% Initial odd H-field, column repeat 
HyoO_tilda=meshgrid(hyoO_tilda,1:Ny).'; 

% Include the reflection coefficient 
Hy0_tilda=0.5*(1+gamma).*Hye0_tilda+0.5*(1-gamma).*HyoO_tilda; 

Hy_tilda=HyO_tilda.*H3D; % Hy_tilda(0+,ky,kz) at x=0 
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Exp2=exp(i*kx*Dx) ; 


Marching range 


o, 

o 


o 

% 3D Parabolic Equation Algorithm 

^ifififififififi!:ifififififififi!:ififififif-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 

o 

C=0; % Setting the range counter 

c=0; % Setting the index counter 


for x=0:Dx:D-Dx; 
c=c+l; 


% Define the total steps based on Dx and D 
% Starting the index counter 


Hy_tilda_Dx=Hy_tilda.*Exp2; % Hy_tilda(Dx,ky,kz) 

% ifft2 wrt to ky and kz 
Hy=Dkz*Dky*(Ny*Nz*ifft2(Hy_tilda_Dx))/(2*pi)^2; 


Hy_H=Hy.*H3D; % Apply the Hamming window in spatial domain 

9''k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

% Insert Screens (Buildings) 

S'if'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

C=C+Dx; % Range Counter 

while M(M==C) % If location of screen = marching distance 

Hy_H(1:round((Hk(c)+Hbs)/Dz),1:round(Wk/Dy/2))=0; 

Hy_H(1:round((Hk(c)+Hbs)/Dz),(Ny-round(Wk/Dy/2)):Ny)=0; 
break 

end 

o 

% Field Strength Convert from 3D to 2D 

9'ifififififififififififififif'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'kif'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

% Measure the field strength at the top of each screen 

9''k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 

J=4; % First column final H-field 

rho=sqrt(C^2+(J)^2); % Distance from transmitter to the receiver 
R=sqrt(((Hk(c)+Hbs)-Ht)^2 + rho^2); 

% Distance from transmitter to receiver 

theta=asin(rho/R) ; % Angle Theta measured relative to the z-axis 

phi=asin(J/rho); % Angle Phi measured relative to the x-axis 

% Free-space g_tilda 

g_tilda_fs=exp(-((cos(theta)*ko).^2*sigmaz^2)/2); 

% Free-space H-field 

Hyfs=-i*ko*sin(theta)*cos(phi).*g_tilda_fs.*exp(i*ko*R)./(4*pi*R); 

% z > 0 of the first column of the final H-field 
A=Hy_H(round( (Hk(c)+Hbs)/Dz)+1,4) ; 

Frel(c)=20*logl0(abs(A/Hyfs)); % Relative Prop Factor 

% Convert propagation factor (F) to E-field by multiplying it by 
% 1/sqrt(ko*R) 

Field(c)=20*logl0(abs(A/Hyfs)/sqrt(ko*R)); % E-Field strength 

FFreeSpace(c)=20*logl0(1/sqrt(ko*R)); % 2D free-space pathless 
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o 


Hyz1=Hy_H(l:Nz/2 + l, :) ; 

HyzO = flipud(Hy_H(2:Nz/2, :) ) ; 

Hye=[Hyzl; HyzO]; 

Hyo=[Hyzl; -HyzO]; 

Hye_tilda=(Dz*Dy)*fft2(Hye) ; 
Hye_tilda_H=Hye_tilda.*H3D; 

Hyo_tilda=(Dz*Dy)*fft2(Hyo) ; 
Hyo_tilda_H=Hyo_tilda.*H3D; 

% Apply reflection 
Hye_tilda_g=0.5*(1+gamma).*Hye_ti 


H-field for z > 0 
H-field for z < 0 
Even H-field 
Odd H-field 

Take the EFT of the even H-tilda 
the Hanning window in freq domain 

Take the EFT of the odd H-tilda 
The Hanning window in freq domain 
oefficient to the even field 
.a_H ; 


% Apply reflection coefficiet to the odd field 
Hyo_tilda_g=0.5*(1-gamma).*Hyo_tilda_H; 


Hy_tildal=Hye_tilda_g + Hyo_tilda_g; % Hy_tilda(x,ky,kz) 


Hy_tilda=Hy_tildal.*H3D; % The Hanning window in freq domain 


end 

save H:\THESIS_2001\MATLAB_CODES\FieldStl00 Field -ASCII 
save H:\THESIS_2001\MATLAB_CODES\Free_Space FFreeSpace -ASCII 

o 

% Plot Results 

o 

figure (1) 

stem(xs2,Hk+Hbs),grid 
xlabel (' Screen Placement (m) ') 
ylabel( 'Screen Height (m)' ) 

figure(2) 

Z=(-1000:Dns:Dist+Dns) 

plot(Z,Frel),grid 

axis([-1000 5000 -80 10]) 

ylabel (' Relative Field Strength (dB) ') 

xlabel (' Screen Position (m) ') 

title([ ' Wk = ' num2str(Wk) ' meters' ' Dx = 'num2str(Dx) ' meters' ' 
Ht = ' num2str(Ht) ' meters']) 

figure(3) 

Z=(-1000:Dns:Dist+Dns)'; 

save H:\THESIS_2001\MATLAB_CODES\RHPlacement Z -ASCII 

plot(Z,Field,Z,FFreeSpace),grid 

axis([-1000 5000 -120 -30]) 

ylabel (' Field Strength (dB) ') 

xlabel (' Screen Position (m) ') 

title([ ' Wk = ' num2str(Wk) ' meters' ' Dx = 'num2str(Dx) ' meters']) 

% Stop stopwatch 
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toe 
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APPENDIX E: A 3D PE MODEL FOR TWO HILLS OF 

SINUSOIDAL SHAPE 


Below is the MatLab program listing of the 3D parabolic equation model use to 
compute the field strengths for two sinusoidal hills with multiple absorbing screens of 
uniform heights and equal spacing. These screens represent a row of buildings or houses 
in residential areas of a well built-up city. The operating frequency is 900 MHz. 

o 

% Declaring Contants 

o 


hold off; 

tic 

c=3e8; 

f=0.9e9; 

Lamda=c/f; 

ko=(2*pi)/Lamda; 

S=400; 

Er=10; 


% Speed of light 
% Operating frequency 
% Wavelength 
% Wavenumber 

% Ground conductivity in mS/m 
% Relative dielectric constant 


Dz=Lamda; 

Dy=Lamda; 

sigmaz=Lamda; 


Delta Z 
Delta Y 
Current 


in meter 
in meter 

source standard deviation 


o 

% Input Parameters 

o 


Ht=57; 
Hbs=7; 


% Transmitting antenna height 
% Receiving antenna height above KE's 


Wk=100; % 

Dns=50; % 

xs=-1500:Dns:4500+Dns; % 

Hk=Hbs+25+25*cos(2*pi*xs/3000) ; 


Knife edge width of 25 meter 
Screen separation 

Horizontal hill width 

% Equal space screen height 


NS=length(xs); 

M=Dns:Dns:Dns*(NS) ; 


% Number of screens 
% Location of the screens 


D=Dns*(NS); 
Dx=5 0 ; 


% Range from the transmitter 
% Incremental range (Dx) in meter 


Ny=1024; 

Nz=1024; 


Sample size 
Sample size 


in the y-direction 
in the z-direction 


^ifififif-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 

o 

% Define range in y-direction and z-direction 


% maximum range in y-direction 


77 


ymax=Ny*Dy/2; 


yl=0:Dy:ymax; 

y=[yl zeros(l,Ny/2-l)]; 


% First half of y range 
% y array 


zmax = Nz*Dz/2; 

zl=0:Dz:zmax; 

z=[zl zeros(l,Nz/2-l)]; 


% Find maximum z range (Zmax) 

% First half of the range of Zmax+1 
% Construct a full z array 


o 

% Define Wavenumbers in Spatial Domain 

^■k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 


kymax=pi/Dy; 

% Maximum wavenumber 

in y-direction 

Dky=2 *kymax/Ny; 

% detla ky 


kyl=-kymax:Dky:kymax; 

% Range of ky 


ky=[kyl(:,(Ny/2)+1:Ny+1) 

kyl(:,2:Ny/2)]; 


kzmax=pi/Dz; 

% Maximum wavenumber 

in z-direction 

Dkz=2 *kzmax/Nz; 

% detla kz 


kzl=-kzmax:Dkz:kzmax; 

% Range of kz 


kz=[kzl(:,(Nz/2)+1:Nz+l) 

kzl ( :,2:Nz/2)]; 



5~kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

o 

% Compute Wavenumbers in Frequency Domain, kx 

^kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 


Ky=meshgrid(ky,1:Nz) ; 
Kz=meshgrid(kz,1:Ny) . ' ; 
repeat 

kx=sqrt (ko''2-Ky. ''2-Kz . ''2) ; 

clear Ky; 
clear Kz; 


Create a NyxNz matrix 
Create a NyxNz matrix 


of ky row-repeat 
of kz column- 


Compute theNy x Nz kx matrix 


Clear matrices Ky and Kz from memory 
To make the algorithm run faster 


Q-kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

o 

% Compute reflection coefficient 

^kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

o 

Erc=Er+i*18*S/(f/le6); % Complex dielectric 

Zs=l/sqrt(Ere); % Impedance 

Gamma=(kz-ko*Zs)./(kz+ko*Zs); % Complex reflection coefficient 

%gamma=-ones(Nz,Ny); % Setting gamma = -1 

gamma=meshgrid(Gamma,1:Ny); % Create reflection coefficient matrix 

% by column repeating 

S'kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

o 

% Define the Hanning Window 

^-kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

o 

Hya=[]; % An empty vector 

for t=0:Ny/2; % Define Ny/2 + 1 points 

if (t>=0 & t<=3*Ny/8) % Define first 3Ny/8 +1 points 

hy=l; 

elseif (t>=3*Ny/8 & t<=Ny/2)% Define next Ny/8 points 
hy=(sin(4*pi*t/Ny))^2; 

end 
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Hya=[Hya hy]; 

end 


% Construct the first half of Hanning 
% window of l+Ny/2 elements 


Yy=fliplr(Hya(:,2:Ny/2) ) ; 
second element to 
HY=[Hya Yy]; 


% Flip left to right of the Hanning 
% windows from the 

% Ny/2 element for total of [lx(Ny/2-l)] 
% Hanning window in y-direction 


Hzb= [ ]; % 

for t=0:Nz/2; % 

if (t>=0 & t<=3*Nz/8) % 

h z = 1 ; 

elseif (t>=3*Nz/8 & t<=Nz/2)% 
hz=(sin(4*pi*t/Nz))^2; 

end 

Hzb=[Hzb hz]; % 

end % 


An empty vector 
Define Nz/2 + 1 points 
Define first 3Nz/8 +1 points 

Define next Nz/8 points 

Construct the first half of Hanning 
window of l+Nz/2 elements 


Yz=fliplr(Hzb(:,2:Nz/2)) ; 


HZ=[Hzb Yz]'; 


% Flip left to right of the Hanning 
% windows from the second element to 
% Nz/2 element for total of [lx(Nz/2-l)] 
% Hanning window in z-direction 


Hmy=meshgrid(HY,1:Nz) ; 
Hmz=meshgrid(HZ,1:Ny) ' ; 


% Row repeat 
% Column repeat 


H3D=(Hmy.*Hmz); 


% 3D Hanning window 


o 

% Gaussian Current Source 

^-'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

o 


f = (1/ (sqrt (2*pi) *sigmaz) ) *exp (- (z-Ht) .'' 2 / (2*sigmaz^2) ) ; 
f_tilda=Dz*fft(f); % Fourier transform of f(z) 

o 

% Initial H-field at x=0 

^■k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 

o 


g_tilda=exp(-kz.^2*sigmaz^2/2); % initial g tilda 

hyeO_tilda=g_tilda.*cos(kz*Ht); % Initial hye_tilda(Of,ky,kz) field 

hyoO_tilda=-i*g_tilda.*sin(kz*Ht); % Initial hyo_tilda(Of,ky,kz) field 

% Initial even H-field, column repeat 
HyeO_tilda=meshgrid(hyeO_tilda,1:Ny).'; 

% Initial odd H-field, column repeat 
HyoO_tilda=meshgrid(hyoO_tilda,1:Ny).'; 

% Include the reflection coefficient 
Hy0_tilda=0.5*(1+gamma).*Hye0_tilda+0.5*(1-gamma).*HyoO_tilda; 

Hy_tilda=HyO_tilda.*H3D; % Hy_tilda(Of,ky,kz) at x=0 

Exp2=exp(i*kx*Dx); % Marching range 
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o 

% 3D PE Basic Algorithm 

^ifififififififififififififififififififififif-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 

o 


C=0; 
c=0 ; 


% Setting the range counter 
% Setting the index counter 


for x=0:Dx:D-Dx; 
D 

c=c+l; 


% Define the total steps based on Dx and 
% Starting the index counter 


Hy_tilda_Dx=Hy_tilda.*Exp2; % Hy_tilda(Dx,ky,kz) 

% ifft2 wrt to ky and kz 
Hy=Dkz*Dky*(Ny*Nz*ifft2(Hy_tilda_Dx))/(2*pi)^2; 


Hy_H=Hy.*H3D; 


% Apply the Hamming window in spatial 


domain 


^ifififif-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 

o 

% Insert Screens (Buildings) 

^■k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k-k 


C=C+Dx; 


% Range Counter 


while M(M==C) % If location of screen = marching distance 

% only work if Dx = Dns; otherwise rewrite 
Hy_H(1:round((Hk(c))/Dz),l:round(Wk/Dy/2))=0; 

Hy_H(1:round((Hk(c))/Dz),(Ny-round(Wk/Dy/2)):Ny)=0; 
break 

end 

^kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

o 

% Field Strength Convert from 3D to 2D 

^kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 


J=4; 


% First column final H-field 
% Distance from Tx to Base Rx 


rho=sqrt(C^2+(J)^2); 

R=sqrt(((Hk(c))-Ht)^2 + rho^2); % Distance from Tx to Rx 

theta=asin(rho/R); % Angle Theta 

phi=asin(J/rho); % Angle Phi 

% Free-space g_tilda 

g_tilda_fs=exp (- ( (cos (theta) *ko) . ''2*sigmaz^2) /2) ; 

% Free-space H-field 

Hyfs=-i*ko*sin(theta)*cos(phi) .*g_tilda_fs.*exp(i*ko*R) ./(4*pi*R) ; 

% z > 0 of the first column of the final H-field 
A=Hy_H(round( (Hk(c) )/Dz)+1,4) ; 


Frel(c)=20*logl0(abs(A/Hyfs)) ; 

Field(c)=2 0*logl0(abs(A/Hyfs)/sqrt(ko*R)) ; 

FFreeSpace(c)=20*logl0(abs(Hyfs)/sqrt(ko*R) ) ; 


% Field strength 
% Field strength 


% H-field for z > 0 


Hyz1=Hy_H(l:Nz/2 + l, :) ; 

HyzO = flipud(Hy_H(2:Nz/2, :) ) ; 


% H-field for z < 0 

% Even H-field 
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Hye=[Hyzl; HyzO] 


Hyo=[Hyzl; -HyzO]; 

Hye_tilda=(Dz*Dy)*fft2(Hye) ; 
Hye_tilda_H=Hye_tilda.*H3D; 

Hyo_tilda=(Dz*Dy)*fft2(Hyo) ; 
Hyo_tilda_H=Hyo_tilda.*H3D; 

% Apply reflect! 
Hye_tilda_g=0.5*(1+gamma).*Hye_ 

% Apply reflect! 
Hyo_tilda_g=0.5*(1-gamma).*Hyo_ 


% Odd H-field 

% Take the FFT of the even H-tilda 
% the Hanning window in freq domain 

% Take the FFT of the odd H-tilda 
% The Hanning window in freq domain 
n coefficient to the even field 
ilda_H; 

n coefficiet to the odd field 
ilda_H; 


Hy_tildal=Hye_tilda_g + Hyo_tilda_g; % Hy_tilda(x,ky,kz) 


Hy_tilda=Hy_tildal.*H3D; % The Hanning window in freq domain 


end 

save H:\THESIS_2001\Base_Line_Codes\TwoHFieldl00 Field -ASCII 


o 

% Plot Results 

^kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

o 

figure (1) 
stem(xs,Hk),grid 
axis([-2000 5000 0 60]) 
xlabel (' Screen Placement (m) ') 
ylabel( 'Screen Height (m)' ) 

figure(2) 

Z=(-1500:Dns:4500+Dns) 

plot(Z,Frel),grid 

axis([-2000 5000 -80 10]) 

ylabel (' Relative Field Strength (dB) ') 

xlabel (' Screen Position (m) ') 

title([ ' Wk = ' num2str(Wk) ' meters' ' Dx = 'num2str(Dx) ' meters' ' 
Ht = ' num2str(Ht) ' meters']) 

figure(3) 

Z=(-1500:Dns:4500+Dns)'; 
plot(Z,Field),grid 
axis([-2000 5000 -120 -30]) 
ylabel (' Field Strength (dB) ') 
xlabel (' Screen Position (m) ') 

title([ ' Wk = ' num2str(Wk) ' meters' ' Dx = 'num2str(Dx) ' meters']) 


%iegend('Parabolic','KnifedB') 

%saveas(gcf, '3D_PE_4RAY.fig') ; 

%saveas(gcf,'3D_PE_4RAY.bmp'); 

toe % Stop stopwatch 
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